1. Description

This R Markdown document describes the analyses performed for the manuscript entitled “Environmental pollution correlates with parasite infection across a riverine landscape” by Io S. Deflem, Seppe Marchand, Federico C.F. Calboli, Joost A.M. Raeymaekers, Filip A.M. Volckaert and Pascal I. Hablützel.

The analyses were run in R 4.2.2

2. Study area and sampling

Up to thirty 0+ three-spined sticklebacks were sampled at 37 locations in the Dijle and Demer basins in Belgium during autumn 2016 under a permit of the Flemish Agency Nature and Forest. Both basins together cover a continuous surface area of 3,624 km² with the furthest two sampling sites being located 117 km apart (distance measured along rivers). All locations included small and relatively slow flowing streams (drop off from highest to lowest point is 18 m) covering a wide range of ecological, hydromorphological, and physico-chemical characteristics. Fish were caught using a dip net.

3. Setting up working environment

# Empty environment
rm(list=ls())

# Print version of R
cat("This script was run with:", version[['version.string']], "\n")
## This script was run with: R version 4.1.2 (2021-11-01)
# Set working directory to location where script is stored
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) # requires installation of package "rstudioapi"

# Loading required libraries
require(BAS)
cat("BAS version", getNamespaceVersion("BAS"), "\n")
## BAS version 1.6.4
require(boral)
cat("boral version", getNamespaceVersion("boral"), "\n")
## boral version 2.0
require(car)
cat("car version", getNamespaceVersion("car"), "\n")
## car version 3.1-2
require(corrplot)
cat("corrplot version", getNamespaceVersion("corrplot"), "\n")
## corrplot version 0.92
require(ggplot2)
cat("ggplot2 version", getNamespaceVersion("ggplot2"), "\n")
## ggplot2 version 3.4.2
require(gplots)
cat("gplots version", getNamespaceVersion("gplots"), "\n")
## gplots version 3.1.3
require(vegan)
cat("vegan version", getNamespaceVersion("vegan"), "\n")
## vegan version 2.6-4
require(factoextra)
cat("factoextra version", getNamespaceVersion("factoextra"), "\n")
## factoextra version 1.0.7.999
require(piecewiseSEM)
cat("piecewiseSEM version", getNamespaceVersion("piecewiseSEM"), "\n")
## piecewiseSEM version 2.1.0
require(remotes)
cat("remotes version", getNamespaceVersion("remotes"), "\n")
## remotes version 2.4.2
require(nlme)
cat("nlme version", getNamespaceVersion("nlme"), "\n")
## nlme version 3.1-155

4 Loading and preparing host and parasite data

Fish were euthanized with a lethal dose of MS222 on the day of capture, following directions of the KU Leuven Animal Ethics Commission, and stored at -20 °C. Fish were kept in separate containers per site at all times. Lab based parasite screening of thawed fish involved placing individual fish in 5 or 10 ml cryo-tubes with 1 or 2 ml of distilled water. Following a vigorous shake of 10 s, the liquid was poured into a Petri dish and ectoparasites were identified and counted using a stereomicroscope. Fish were rinsed and checked again for the presence of ectoparasites on skin and fins. The intestines were examined for endoparasites. Before dissection, fish weight (± 0.1 mg) and standard length (± 1 mm) were recorded. To quantify body condition, we calculated the scaled mass index (SMI; Maceda-Veiga et al., 2014; Peig & Green, 2009). Sex was determined during dissection by inspection of gonad development. A total of 668 fish were dissected, which amounts to approximately 20 fish per location, with the exception of seven locations where only 10 fish were screened for the presence of macroparasites. Ecto- and endoparasites were morphologically identified to species level whenever possible.

# Parasite data
data <- read.csv("data_2016_2303.csv", sep=';')
data$site <- as.factor(data$site)

# Calculate parasite parameters
#names(data)
# Parasite data is overdispersed (mostly so for Trichodina), if using average abundance data, species matrix needs to be transformed
# Remove individual fish for which no parasites were counted
datao <- na.omit(data[,c(1,22:24,26:32)]) 

ddata <- dispweight(datao[,-1]) # Correct for overdispersion of the parasite count data
avab <- aggregate(ddata, by = list(datao[, "site"]), function(x){mean(x, na.rm =T)}) # Calculate average abundances per site
prev = aggregate(data[,c(22:24,26:32)], by = list(data[, "site"]), function(x){sum(x >0, na.rm = T)/length(x)}) # Calculate prevalence per site
medin = aggregate(data[,c(22:24,26:32)], by = list(data[, "site"]), function(x){median(x[x >0], na.rm = T)}) # Calculate median infection intensity per site
pa = aggregate(data[,c(22:24,26:32)], by = list(data[, "site"]), function(x){ifelse(mean(x, na.rm =T)>0,1,0)}) # Assess presence/absence per site

avab[is.na(avab)] <- 0
prev[is.na(prev)] <- 0
medin[is.na(medin)] <- 0

# Host condition data
avcondition <- aggregate(data$SMI, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2] # Calculate average host condition (SMI) per site
avlength <- aggregate(data$length, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2] # Calculate average host length per site

# Parasite index
sgyr <- 1:nrow(data)
stri <- 1:nrow(data)
sglu <- 1:nrow(data)
scon <- 1:nrow(data)
scysl <- 1:nrow(data)
spro <- 1:nrow(data)
saca <- 1:nrow(data)
scam <- 1:nrow(data)
sang <- 1:nrow(data)

for(j in 1:nrow(data)){
  sgyr[j] <- data$Gyr[j]/sd(data$Gyr, na.rm=T)
  stri[j] <- data$Tri[j]/sd(data$Tri, na.rm=T)
  sglu[j] <- data$Glu[j]/sd(data$Glu, na.rm=T)
  scon[j] <- data$Con[j]/sd(data$Con, na.rm=T)
  scysl[j] <- data$CysL[j]/sd(data$CysL, na.rm=T)
  spro[j] <- data$Pro[j]/sd(data$Pro, na.rm=T)
  saca[j] <- data$Aca[j]/sd(data$Aca, na.rm=T)
  scam[j] <- data$Cam[j]/sd(data$Cam, na.rm=T)
  sang[j] <- data$Ang[j]/sd(data$Ang, na.rm=T)
}

PI <- 1:nrow(data)
for(j in 1:nrow(data)){
  PI[j] <- 10/max(sgyr, na.rm=T)*sgyr[j] + 10/max(stri, na.rm=T)*stri[j] + 10/max(sglu, na.rm=T)*sglu[j] +
    10/max(scon, na.rm=T)*scon[j] + 10/max(scysl, na.rm=T)*scysl[j] + 10/max(spro, na.rm=T)*spro[j] +
    10/max(saca, na.rm=T)*saca[j] + 10/max(scam, na.rm=T)*scam[j] + 10/max(sang, na.rm=T)*sang[j]
}

PI_ecto <- 1:nrow(data)
for(j in 1:nrow(data)){
    PI_ecto[j] <- 10/max(sgyr, na.rm=T)*sgyr[j] + 10/max(stri, na.rm=T)*stri[j] + 10/max(sglu, na.rm=T)*sglu[j]
}

PI_endo <- 1:nrow(data)
for(j in 1:nrow(data)){
  PI_endo[j] <-     10/max(scon, na.rm=T)*scon[j] + 10/max(scysl, na.rm=T)*scysl[j] + 10/max(spro, na.rm=T)*spro[j] +
    10/max(saca, na.rm=T)*saca[j] + 10/max(scam, na.rm=T)*scam[j] + 10/max(sang, na.rm=T)*sang[j]
}

avPI <- aggregate(PI, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2] # Calculate average parasite index per site
avPI_ecto <- aggregate(PI_ecto, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2] # Calculate average parasite index for ectoparasites per site
avPI_endo <- aggregate(PI_endo, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2] # Calculate average parasite index for endoparasites per site

write.csv(prev, "supplementary_table_2.csv", row.names=FALSE)

5 Loading and preparing environmental and spatial data

Physico-chemical data was provided by the Flemish Environmental Agency (VMM). Each fish sampling site was chosen at or near an environmental monitoring site of VMM. Water parameters include water temperature, pH, conductivity, dissolved oxygen (O2), saturation with dissolved oxygen, and Biochemical and Chemical Oxygen Demand (BOD and COD). Nutrient related water parameters include measurements of nitrate (NO3-), nitrite (NO2-), Kjeldahl nitrogen (KjN), total nitrogen (Nt), ammonium (NH4+), and total phosphorus (Pt). Following removal of strong collinear variables (significant correlation with P < 0.05 and Pearson correlation coefficient > |0.6|; Dormann et al., 2013), six environmental physico-chemical variables were retained (temperature, conductivity, COD, saturation with dissolved oxygen, ammonium, and total nitrogen), representing different aspects of water quality. For each water parameter, the average value of the year before sampling was calculated based on monthly monitoring data. Additionally, two hydromorphological variables were included: Tthe presence of a pool-riffle pattern and meanders were noted during field sampling and these parameters were included as binary variables (presence/absence) for a representation of abiotic habitat structure. Spatial (waterway) distances were calculated using the Network Analyst toolbox in ArcGIS. Upstream distance was defined as the maximal upstream distance from the sampling location. Network peripherality was calculated as the average waterway distance of a sampling location to all other locations. Hence, a total of eight environmental and two spatial variables were included in the statistical analysis.

# Environmental data (VMM)
environment <- read.csv("Environment_update.csv", sep=';')
# Spatial variables: network centrality and upstream distance
spavar <- read.csv("space2.csv", sep=';')
#plot(spavar$netcen); plot(density(spavar$netcen))
#plot(spavar$updist); plot(density(spavar$updist))

# Environmental data (from field observations)
field_data <- read.csv("field_data.csv", sep=',')
environment2 <- cbind(environment[,c(1,49,52:53,55,57,60,63)], field_data[-c(8,10,25,27,37),33:34], spavar[,c(2,3)])
environment2$pool_riffle <- as.factor(environment2$pool_riffle)
environment2$meander <- as.factor(environment2$meander)

netcen <- spavar$netcen
updist <- spavar$updist

supplementary_table_1 = cbind(environment2, netcen, updist)
write.csv(supplementary_table_1, "supplementary_table_1.csv", row.names=FALSE)

We used univariate generalized linear models to investigate how landscape-level effects modify infection patterns of individual parasite taxa, host size and condition. We kept the statistical models linear (as opposed to polynomial) and only considered main effects (i.e. no interaction terms) because we had no prior information from this study system that more complex models were necessary and because the study design with (only) 37 sampling sites was not intended for non-linear models. Ten explanatory variables (temperature, conductivity, COD, saturation with dissolved oxygen, ammonium, total nitrogen, the presence of pool-riffle patterns and meanders, upstream distance, and network peripherality) were included.

6. Univeriate analysis using Bayesian Model Averaging

Univariate analyses - We used generalized linear models in a BMA approach to understand how infection with individual parasite taxa relate to host characteristics (length and condition), environmental conditions as well as spatial distance among sampling sites. Parasite infection was calculated in three ways at the host population level: average abundance (mean parasites per host), prevalence (percentage of infected hosts) and median infection intensity (median number of parasites in infected hosts). We calculated the individual parasitation index (IPI) following Kalbe et al. (2002) as a measurement for total parasite abundance and species richness for each individual fish. This index was calculated for all parasite species combined, and for ecto- and endoparasite species separately. For these models, we assumed a normal error distribution (which appeared to be justified, see Supplementary Figures S1-S2) and applied a Jeffrey-Zellner-Siow prior. Model assumptions (homoscedasticity of the variances and normal distribution of the errors) were assessed using the generic model plot function in R and did onlynot clearly deviate in any of the models for rare parasites. We followed a normal distribution, and not Poisson or negative binomial, for the parasite data for the common species (Trichodina sp. and Gyrodactylus spp.) and the individual parasitation index as the parameters used are deviates from count data. Rare parasites (Glugea, Contracaecum, Anguillicoloides, and unidentified cysts) were excluded from the univariate analysis because there was not enough data to obtain a good fit of the models.For rare parasites (Contracaecum sp. and Anguillicoloides crassus), we used population-level presence-absence data assuming a binomial error distribution and a uniformly distributed BIC prior. Due to low prevalences, the other parasites were not included in the species-specific models. Explanatory variables were considered important when they had a posterior inclusion probability (PIP) of 0.5. To account for overdispersion in the parasite counts, we transformed the data by downweighting overdispersed taxa following Clarke et al. (2006) using the dispweight function in the R package vegan v2.5.6 (Oksanen et al., 2013).

6.0.1 Figure 2

# Make a matrix for PIP (Posterior Inclusion Probability)
PIP <- matrix(nrow=12, ncol=14)
rownames(PIP) <- c("Host condition", "Host length", "Temperature", "Oxygen saturation", "Conductivity", "COD", "Ammonium", "Total nitrogen", "Pool riffle pattern", "Meander", "Network peripherality", "Upstream distance")
colnames(PIP) <- c("Host condition", "Host size", "Gyrodactylus abundance", "Gyrodactylus prevalence", "Gyrodactylus infection intensity", "Trichodina abundance", "Trichodina prevalence", "Trichodina infection intensity", "Glugea", "Contracaecum", "Aguillicola",
                   "PI", "PI ecto", "PI endo")  
#Condition
bas.model <- bas.lm(avcondition ~ T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
plot(bas.model)

coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(3:12),1] <- pip[2:11,1]*sign(coef.model$postmean[2:11])

#Length
bas.model <- bas.lm(avlength ~ T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(3:12),2] <- pip[2:11,1]*sign(coef.model$postmean[2:11])

#Gyrodactylus abundance
bas.model <- bas.lm(avab$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#Gyrodactylus prevalence
bas.model <- bas.lm(prev$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),4] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#Gyrodactylus infection intensity
bas.model <- bas.lm(medin$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),5] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#Trichodina abundance
bas.model <- bas.lm(avab$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),6] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#Trichodina prevalence
bas.model <- bas.lm(prev$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),7] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#Trichodina infection intensity
bas.model <- bas.lm(medin$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),8] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#Glugea
bas.model <- bas.glm(pa$Glu ~  avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, betaprior=g.prior(100), family=binomial)
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),9] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#Contracaecum
bas.model <- bas.glm(pa$Con ~  avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, betaprior=g.prior(100), family=binomial)
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),10] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#Anguillicola
bas.model <- bas.glm(pa$Ang ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, betaprior=g.prior(100), family=binomial)
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),11] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#PI
bas.model <- bas.lm(avPI ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),12] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#PI ecto
bas.model <- bas.lm(avPI_ecto ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),13] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

#PI endo
bas.model <- bas.lm(avPI_endo ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
coef.model <- coef(bas.model)
pip <- summary(bas.model)
PIP[c(1:12),14] <- pip[2:13,1]*sign(coef.model$postmean[2:13])
x = round(PIP, digits=2)
x[abs(PIP)<0.5] <- ""
x[abs(PIP)>0.5] <- "+"
heatmap.2(PIP[,-c(9,10,11)],
          cellnote = x[,-c(9,10,11)],
          #main = "Correlation",
          notecex=1,
          notecol="white",
          density.info="none",
          trace="none",
          margins =c(10,8),
          col=redblue(256),
          dendrogram="both",
          cexRow = 0.7,
          cexCol = 0.7,
          key.title = "PIP",
          lhei = c(1,3),
          lwid = c(0.5, 0.5),
          srtCol = 45,
          labCol = c("Host condition", "Host length",
                     expression(paste(italic("Gyrodactylus"),  " spp. - abundance")),
                     expression(paste(italic("Gyrodactylus"), " spp. - prevalence")),
                     expression(paste(italic("Gyrodactylus"), " spp. - infection intensity")),
                     expression(paste(italic("Trichodina"), " sp. - abundance")),
                     expression(paste(italic("Trichodina"), " sp. - prevalence")), 
                     expression(paste(italic("Trichodina"), " sp. - infection intensity")), 
                     expression("I"[PI]*" all parasites"), 
                     expression("I"[PI]*" ectoparasites"), 
                     expression("I"[PI]*" endoparasites")) ,
          #Colv="NA"
          ) 

pdf("Figure2.pdf", width = 7.29, height = 4.5)
heatmap.2(PIP[,-c(9,10,11)],
          cellnote = x[,-c(9,10,11)],
          #main = "Correlation",
          notecex=1,
          notecol="white",
          density.info="none",
          trace="none",
          margins =c(10,8),
          col=redblue(256),
          dendrogram="both",
          cexRow = 0.7,
          cexCol = 0.7,
          key.title = "PIP",
          lhei = c(1,3),
          lwid = c(0.5, 0.5),
          srtCol = 45,
          labCol = c("Host condition", "Host length", 
                     expression(paste(italic("Gyrodactylus"),  " spp. - abundance")), 
                     expression(paste(italic("Gyrodactylus"), " spp. - prevalence")), 
                     expression(paste(italic("Gyrodactylus"), " spp. - infection intensity")), 
                     expression(paste(italic("Trichodina"), " sp. - abundance")), 
                     expression(paste(italic("Trichodina"), " sp. - prevalence")), 
                     expression(paste(italic("Trichodina"), " sp. - infection intensity")), 
                     expression("I"[PI]*" all parasites"), 
                     expression("I"[PI]*" ectoparasites"), 
                     expression("I"[PI]*" endoparasites")) ,
          #Colv="NA"
          ) 

dev.off()
## png 
##   2

6.1 Variation in host condition

bas.model <- bas.lm(avcondition ~  T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + 
                      pool_riffle + meander + netcen + updist, 
                    data=environment2, prior="JZS")
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y) model 1    model 2    model 3    model 4    model 5
## Intercept       1.00000000  1.0000  1.0000000  1.0000000  1.0000000  1.0000000
## T_av            0.12386945  0.0000  1.0000000  0.0000000  0.0000000  0.0000000
## O2_sat_av       0.05073157  0.0000  0.0000000  0.0000000  0.0000000  0.0000000
## Con_av          0.03919858  0.0000  0.0000000  0.0000000  0.0000000  0.0000000
## COD_av          0.05197281  0.0000  0.0000000  0.0000000  0.0000000  1.0000000
## NH4._av         0.04846220  0.0000  0.0000000  0.0000000  0.0000000  0.0000000
## Nt_av           0.05248386  0.0000  0.0000000  0.0000000  0.0000000  0.0000000
## pool_riffle1    0.06729856  0.0000  0.0000000  0.0000000  1.0000000  0.0000000
## meander1        0.06527705  0.0000  0.0000000  0.0000000  0.0000000  0.0000000
## netcen          0.06898689  0.0000  0.0000000  1.0000000  0.0000000  0.0000000
## updist          0.04496936  0.0000  0.0000000  0.0000000  0.0000000  0.0000000
## BF                      NA  1.0000  0.6511099  0.3462387  0.2957714  0.2212089
## PostProbs               NA  0.6912  0.0450000  0.0239000  0.0204000  0.0153000
## R2                      NA  0.0000  0.0906000  0.0565000  0.0478000  0.0315000
## dim                     NA  1.0000  2.0000000  2.0000000  2.0000000  2.0000000
## logmarg                 NA  0.0000 -0.4290769 -1.0606268 -1.2181684 -1.5086476
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
#abs(coef.model$postmean)-2*coef.model$postsd > 0
plot(confint(coef.model, parm = 2:11))
## Warning in arrows(x[not.deg], ci[not.deg, 1], x[not.deg], ci[not.deg, 2], :
## zero-length arrow is of indeterminate angle and so skipped

## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'condition.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(3:12),1] <- pip[2:11,1]*sign(coef.model$postmean[2:11])
#coef.model$postmean[2:11]

6.2 Variation in host length

bas.model <- bas.lm(avlength ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + 
                      pool_riffle + meander + netcen + updist, 
                    data=environment2, prior="JZS")
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)   model 1    model 2   model 3  model 4   model 5
## Intercept        1.0000000 1.0000000 1.00000000 1.0000000 1.000000 1.0000000
## T_av             0.1914598 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## O2_sat_av        0.1407247 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## Con_av           0.4373524 0.0000000 0.00000000 1.0000000 1.000000 0.0000000
## COD_av           0.1445232 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## NH4._av          0.2026061 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## Nt_av            0.1889236 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## pool_riffle1     0.2259504 0.0000000 0.00000000 0.0000000 0.000000 1.0000000
## meander1         0.3217534 0.0000000 0.00000000 0.0000000 1.000000 0.0000000
## netcen           0.6121912 1.0000000 0.00000000 0.0000000 0.000000 1.0000000
## updist           0.1480063 0.0000000 0.00000000 0.0000000 0.000000 0.0000000
## BF                      NA 0.8278694 0.07277547 0.2894406 1.000000 0.6785609
## PostProbs               NA 0.1582000 0.13910000 0.0553000 0.042500 0.0288000
## R2                      NA 0.2306000 0.00000000 0.1819000 0.314300 0.2982000
## dim                     NA 2.0000000 1.00000000 2.0000000 3.000000 3.0000000
## logmarg                 NA 2.4314765 0.00000000 1.3805712 2.620376 2.2325953
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
#abs(coef.model$postmean)-2*coef.model$postsd > 0
plot(confint(coef.model, parm = 2:11))
## Warning in arrows(x[not.deg], ci[not.deg, 1], x[not.deg], ci[not.deg, 2], :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(x[not.deg], ci[not.deg, 1], x[not.deg], ci[not.deg, 2], :
## zero-length arrow is of indeterminate angle and so skipped

## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'length.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(3:12),1] <- pip[2:11,1]*sign(coef.model$postmean[2:11])
#coef.model$postmean[2:11]

# Prediction plot
newdata = as.data.frame(cbind(rep(mean(environment2$T_av), 37), 
                rep(mean(environment2$O2_sat_av), 37),
                rep(mean(environment2$Con_av), 37),
                rep(mean(environment2$COD_av), 37),
                rep(mean(environment2$NH4._av), 37),
                rep(mean(environment2$Nt_av), 37),
                rep(1, 37),
                rep(1, 37),
                rep(mean(netcen), 37),
                rep(mean(updist), 37)))
colnames(newdata) <- c("T_av", "O2_sat_av", "Con_av", "COD_av", "NH4._av", "Nt_av", "pool_riffle", "meander", "netcen", "updist")
newdata[,"pool_riffle"] <- as.factor(newdata[,"pool_riffle"]); newdata[,"meander"] <- as.factor(newdata[,"meander"])
newdata1 <- newdata; newdata1[,"netcen"] <- netcen
BMA_avlength_netcen <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)

figure3j <- ggplot(environment2, aes(netcen, BMA_avlength_netcen$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_avlength_netcen$fit-BMA_avlength_netcen$se.bma.fit), ymax = (BMA_avlength_netcen$fit+BMA_avlength_netcen$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=netcen, y=avlength)) +
  labs(x=expression("Network peripherality [m]"), y=expression("Average host length [mm]")) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "j)")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
figure3j

6.3 Variation in Gyrodactylus infection

6.3.1 Mean abundance

bas.model <- bas.lm(avab$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)  model 1    model 2   model 3   model 4   model 5
## Intercept        1.0000000 1.000000 1.00000000 1.0000000 1.0000000 1.0000000
## avlength         0.5500651 1.000000 0.00000000 0.0000000 0.0000000 1.0000000
## avcondition      0.3049274 0.000000 0.00000000 0.0000000 0.0000000 1.0000000
## T_av             0.1620852 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## O2_sat_av        0.2017699 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## Con_av           0.2621665 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## COD_av           0.8086300 1.000000 0.00000000 1.0000000 1.0000000 1.0000000
## NH4._av          0.2417809 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## Nt_av            0.8618345 1.000000 1.00000000 1.0000000 1.0000000 1.0000000
## pool_riffle1     0.1760484 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## meander1         0.1775842 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## netcen           0.3006370 0.000000 0.00000000 1.0000000 0.0000000 0.0000000
## updist           0.1648814 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## BF                      NA 1.000000 0.03902505 0.5452757 0.1537678 0.8271829
## PostProbs               NA 0.072100 0.05160000 0.0393000 0.0369000 0.0265000
## R2                      NA 0.524000 0.28790000 0.5059000 0.4100000 0.5631000
## dim                     NA 4.000000 2.00000000 4.0000000 3.0000000 5.0000000
## logmarg                 NA 6.997337 3.75378552 6.3908733 5.1250253 6.8076077
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                       2.5%        97.5%          beta
## Intercept     9.498176e-02 2.044717e-01  1.489134e-01
## avlength     -2.515130e-02 2.264807e-05 -8.122675e-03
## avcondition  -9.805559e-01 0.000000e+00 -1.500326e-01
## T_av         -2.205915e-02 1.634873e-02 -1.493132e-04
## O2_sat_av    -1.678010e-03 5.033937e-03  3.794809e-04
## Con_av       -9.556692e-05 5.266581e-04  5.726233e-05
## COD_av        0.000000e+00 1.300051e-02  6.652983e-03
## NH4._av      -6.709599e-02 7.179572e-02 -1.888173e-04
## Nt_av         0.000000e+00 5.230759e-02  2.930918e-02
## pool_riffle1 -3.796033e-02 1.029158e-01  3.966550e-03
## meander1     -9.234735e-02 5.004225e-02 -4.806259e-03
## netcen       -5.720044e-07 1.018422e-05  1.332021e-06
## updist       -1.325467e-06 1.732882e-06  4.618571e-08
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

6.3.1.1 Prediction plot for marginal effect of host condition on average Gyrodactylus infection

# Prediction plot
newdata = as.data.frame(cbind(rep(mean(avlength), 37),
                rep(mean(avcondition), 37),
                rep(mean(environment2$T_av), 37), 
                rep(mean(environment2$O2_sat_av), 37),
                rep(mean(environment2$Con_av), 37),
                rep(mean(environment2$COD_av), 37),
                rep(mean(environment2$NH4._av), 37),
                rep(mean(environment2$Nt_av), 37),
                rep(1, 37),
                rep(1, 37),
                rep(mean(netcen), 37),
                rep(mean(updist), 37)))
colnames(newdata) <- c("avlength", "avcondition", "T_av", "O2_sat_av", "Con_av", "COD_av", "NH4._av", "Nt_av", "pool_riffle", "meander", "netcen", "updist")
newdata[,"pool_riffle"] <- as.factor(newdata[,"pool_riffle"]); newdata[,"meander"] <- as.factor(newdata[,"meander"])
newdata1 <- newdata; newdata1[,"avcondition"] <- avcondition
BMA_Gyr_avcond <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3k = ggplot(environment2, aes(avcondition, BMA_Gyr_avcond$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_Gyr_avcond$fit-BMA_Gyr_avcond$se.bma.fit), ymax = (BMA_Gyr_avcond$fit+BMA_Gyr_avcond$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=avcondition, y=avab$Gyr)) +
  labs(x=expression("Average host condition"), y=expression(paste(italic("Gyrodactylus"), " spp. - abundance"))) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "k)")

figure3k

6.3.1.2 Prediction plot for marginal effect of COD on average Gyrodactylus infection

newdata1 <- newdata; newdata1[,"COD_av"] <- environment2$COD_av
BMA_Gyr_COD_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3a = ggplot(environment2, aes(COD_av, BMA_Gyr_COD_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_Gyr_COD_av$fit-BMA_Gyr_COD_av$se.bma.fit), ymax = (BMA_Gyr_COD_av$fit+BMA_Gyr_COD_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=COD_av, y=avab$Gyr)) +
  labs(x=expression("Chemical Oxygen Demand [mg/L]"), y=expression(paste(italic("Gyrodactylus"), " spp. - abundance"))) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "a)")
figure3a

6.3.2.3 Prediction plot for marginal effect of total nitrogen on average Gyrodactylus infection

newdata1 <- newdata; newdata1[,"Nt_av"] <- environment2$Nt_av
BMA_Gyr_Nt_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3e = ggplot(environment2, aes(Nt_av, BMA_Gyr_Nt_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_Gyr_Nt_av$fit-BMA_Gyr_Nt_av$se.bma.fit), ymax = (BMA_Gyr_Nt_av$fit+BMA_Gyr_Nt_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=Nt_av, y=avab$Gyr)) +
  labs(x=expression("Nt [mg/L]"), y=expression(paste(italic("Gyrodactylus"), " spp. - abundance"))) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "e)")
figure3e

6.3.2 Median infection intensity

bas.model <- bas.lm(medin$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals

plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y) model 1    model 2   model 3    model 4    model 5
## Intercept       1.00000000  1.0000  1.0000000  1.000000  1.0000000  1.0000000
## avlength        0.02065469  0.0000  0.0000000  0.000000  0.0000000  0.0000000
## avcondition     0.04286905  0.0000  1.0000000  0.000000  0.0000000  0.0000000
## T_av            0.02991678  0.0000  0.0000000  0.000000  0.0000000  0.0000000
## O2_sat_av       0.02400544  0.0000  0.0000000  0.000000  0.0000000  0.0000000
## Con_av          0.02143314  0.0000  0.0000000  0.000000  0.0000000  0.0000000
## COD_av          0.03655302  0.0000  0.0000000  0.000000  1.0000000  0.0000000
## NH4._av         0.03250611  0.0000  0.0000000  0.000000  0.0000000  1.0000000
## Nt_av           0.04088810  0.0000  0.0000000  1.000000  0.0000000  0.0000000
## pool_riffle1    0.02307275  0.0000  0.0000000  0.000000  0.0000000  0.0000000
## meander1        0.02038443  0.0000  0.0000000  0.000000  0.0000000  0.0000000
## netcen          0.02173649  0.0000  0.0000000  0.000000  0.0000000  0.0000000
## updist          0.02738794  0.0000  0.0000000  0.000000  0.0000000  0.0000000
## BF                      NA  1.0000  0.3176354  0.302852  0.2488803  0.2411373
## PostProbs               NA  0.7775  0.0206000  0.019600  0.0161000  0.0156000
## R2                      NA  0.0000  0.0517000  0.049100  0.0381000  0.0363000
## dim                     NA  1.0000  2.0000000  2.000000  2.0000000  2.0000000
## logmarg                 NA  0.0000 -1.1468512 -1.194511 -1.3907833 -1.4223890
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                  2.5%   97.5%          beta
## Intercept    1.564588 3.19689  2.405405e+00
## avlength     0.000000 0.00000 -5.938396e-04
## avcondition  0.000000 0.00000 -2.939823e-01
## T_av         0.000000 0.00000 -7.368692e-03
## O2_sat_av    0.000000 0.00000 -4.581964e-04
## Con_av       0.000000 0.00000 -1.594570e-05
## COD_av       0.000000 0.00000  1.609035e-03
## NH4._av      0.000000 0.00000  1.027436e-02
## Nt_av        0.000000 0.00000  7.570653e-03
## pool_riffle1 0.000000 0.00000  1.051497e-02
## meander1     0.000000 0.00000  2.511911e-03
## netcen       0.000000 0.00000  3.990331e-07
## updist       0.000000 0.00000  4.560731e-07
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

6.3.2 Prevalence

bas.model <- bas.lm(prev$Gyr ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals

plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)  model 1    model 2   model 3   model 4     model 5
## Intercept       1.00000000 1.000000 1.00000000 1.0000000 1.0000000  1.00000000
## avlength        0.16685753 0.000000 0.00000000 0.0000000 1.0000000  0.00000000
## avcondition     0.07728991 0.000000 0.00000000 0.0000000 0.0000000  0.00000000
## T_av            0.07242807 0.000000 0.00000000 0.0000000 0.0000000  0.00000000
## O2_sat_av       0.09667222 0.000000 0.00000000 0.0000000 0.0000000  0.00000000
## Con_av          0.58076129 1.000000 0.00000000 1.0000000 0.0000000  0.00000000
## COD_av          0.17268847 0.000000 0.00000000 1.0000000 0.0000000  0.00000000
## NH4._av         0.09134534 0.000000 0.00000000 0.0000000 0.0000000  0.00000000
## Nt_av           0.11056039 0.000000 0.00000000 0.0000000 0.0000000  0.00000000
## pool_riffle1    0.13922799 0.000000 0.00000000 0.0000000 0.0000000  0.00000000
## meander1        0.09905733 0.000000 0.00000000 0.0000000 0.0000000  0.00000000
## netcen          0.11864979 0.000000 0.00000000 0.0000000 0.0000000  1.00000000
## updist          0.06837623 0.000000 0.00000000 0.0000000 0.0000000  0.00000000
## BF                      NA 1.000000 0.08283414 0.8857497 0.1246735  0.07409002
## PostProbs               NA 0.228800 0.22750000 0.0369000 0.0285000  0.01700000
## R2                      NA 0.233300 0.00000000 0.3039000 0.1341000  0.10730000
## dim                     NA 2.000000 1.00000000 3.0000000 2.0000000  2.00000000
## logmarg                 NA 2.490915 0.00000000 2.3695941 0.4088580 -0.11155941
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                       2.5%        97.5%          beta
## Intercept     2.498255e-01 3.931726e-01  3.224082e-01
## avlength     -2.252212e-02 2.446193e-04 -2.431148e-03
## avcondition  -4.252456e-01 6.599841e-03 -1.986363e-02
## T_av         -1.220894e-02 6.769487e-03 -7.123130e-04
## O2_sat_av    -4.258310e-03 6.559308e-05 -2.662986e-04
## Con_av        0.000000e+00 8.103246e-04  3.112530e-04
## COD_av        0.000000e+00 7.914329e-03  9.376253e-04
## NH4._av      -1.771781e-03 4.657468e-02  1.064409e-03
## Nt_av         0.000000e+00 2.507558e-02  1.612509e-03
## pool_riffle1  0.000000e+00 1.573096e-01  1.465705e-02
## meander1     -1.017710e-01 1.415960e-02 -6.966751e-03
## netcen       -1.412903e-07 7.995999e-06  5.765326e-07
## updist       -7.999739e-07 6.404386e-07 -6.922933e-09
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

6.3.2.3 Prediction plot for marginal effect of conductivity on average Gyrodactylus infection

newdata1 <- newdata; newdata1[,"Con_av"] <- environment2$Con_av
BMA_Gyr_Con_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3i = ggplot(environment2, aes(Con_av, BMA_Gyr_Con_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_Gyr_Con_av$fit-BMA_Gyr_Con_av$se.bma.fit), ymax = (BMA_Gyr_Con_av$fit+BMA_Gyr_Con_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=Con_av, y=prev$Gyr)) +
  labs(x=expression(paste("Conductivity [", mu, "S/cm]")), y=expression(paste(italic("Gyrodactylus"), " spp. - prevalence"))) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "i)")
figure3i

6.4 Variation in Trichodina infection

6.4.1 Mean abundance

bas.model <- bas.lm(avab$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)   model 1  model 2  model 3   model 4   model 5
## Intercept       1.00000000 1.0000000 1.000000 1.000000 1.0000000 1.0000000
## avlength        0.12444042 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## avcondition     0.09559152 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## T_av            0.08500185 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## O2_sat_av       0.09751169 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## Con_av          0.42180029 0.0000000 1.000000 0.000000 0.0000000 0.0000000
## COD_av          0.46724539 0.0000000 1.000000 0.000000 0.0000000 1.0000000
## NH4._av         0.21465122 0.0000000 0.000000 0.000000 1.0000000 0.0000000
## Nt_av           0.34576751 0.0000000 0.000000 1.000000 0.0000000 1.0000000
## pool_riffle1    0.13743699 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## meander1        0.08968689 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## netcen          0.10313787 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## updist          0.08604777 0.0000000 0.000000 0.000000 0.0000000 0.0000000
## BF                      NA 0.0239404 1.000000 0.170643 0.1105942 0.2711811
## PostProbs               NA 0.1650000 0.104400 0.098000 0.0635000 0.0283000
## R2                      NA 0.0000000 0.358500 0.209300 0.1890000 0.3063000
## dim                     NA 1.0000000 3.000000 2.000000 2.0000000 3.0000000
## logmarg                 NA 0.0000000 3.732188 1.964006 1.5303004 2.4272192
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                       2.5%        97.5%          beta
## Intercept     1.266097e-01 2.732195e-01  2.002413e-01
## avlength     -1.608394e-02 3.767937e-04 -1.114390e-03
## avcondition  -4.632070e-01 1.274208e-01 -2.394139e-02
## T_av         -1.438322e-02 1.013193e-02 -1.479738e-04
## O2_sat_av    -1.905880e-03 3.300222e-03  1.132458e-04
## Con_av        0.000000e+00 7.774852e-04  2.109076e-04
## COD_av        0.000000e+00 1.393205e-02  4.203136e-03
## NH4._av      -2.205558e-03 1.033887e-01  1.096330e-02
## Nt_av        -4.771345e-05 4.961465e-02  1.044647e-02
## pool_riffle1  0.000000e+00 1.466687e-01  1.212839e-02
## meander1     -9.597451e-02 4.095486e-03 -1.678368e-03
## netcen       -6.164040e-06 1.239477e-06 -2.846749e-07
## updist       -1.511702e-06 6.054639e-07 -1.915234e-08
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

6.4.1 Median infection intensity

bas.model <- bas.lm(medin$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)  model 1    model 2   model 3   model 4   model 5
## Intercept        1.0000000 1.000000 1.00000000 1.0000000 1.0000000 1.0000000
## avlength         0.2881774 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## avcondition      0.2224163 0.000000 0.00000000 0.0000000 1.0000000 0.0000000
## T_av             0.1574615 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## O2_sat_av        0.1621652 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## Con_av           0.7912882 1.000000 0.00000000 1.0000000 1.0000000 1.0000000
## COD_av           0.7219333 1.000000 0.00000000 0.0000000 1.0000000 1.0000000
## NH4._av          0.3248706 0.000000 1.00000000 1.0000000 0.0000000 0.0000000
## Nt_av            0.3218716 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## pool_riffle1     0.2081402 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## meander1         0.2655250 0.000000 0.00000000 0.0000000 0.0000000 1.0000000
## netcen           0.2107764 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## updist           0.1536937 0.000000 0.00000000 0.0000000 0.0000000 0.0000000
## BF                      NA 1.000000 0.04957628 0.2279222 0.5092549 0.4837322
## PostProbs               NA 0.148800 0.04060000 0.0339000 0.0227000 0.0216000
## R2                      NA 0.449900 0.26810000 0.3987000 0.4816000 0.4799000
## dim                     NA 3.000000 2.00000000 3.0000000 4.0000000 4.0000000
## logmarg                 NA 6.288599 3.28435583 4.8098478 5.6137920 5.5623749
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                       2.5%        97.5%          beta
## Intercept     1.660851e+01 3.395723e+01  2.547297e+01
## avlength     -3.541722e+00 2.164229e-02 -5.158998e-01
## avcondition  -1.441975e+02 4.413419e+00 -1.375058e+01
## T_av         -1.880970e+00 4.198759e+00  2.065105e-01
## O2_sat_av    -5.817500e-01 4.198218e-01 -1.498026e-02
## Con_av        0.000000e+00 1.277170e-01  6.585727e-02
## COD_av        0.000000e+00 1.980322e+00  9.195418e-01
## NH4._av      -1.801730e+00 1.628155e+01  2.222211e+00
## Nt_av        -1.003826e-01 6.309924e+00  1.012166e+00
## pool_riffle1 -4.799011e-01 2.544697e+01  2.247962e+00
## meander1     -2.640241e+01 1.529459e-01 -3.523477e+00
## netcen       -1.457582e-03 1.774935e-04 -1.366046e-04
## updist       -2.281233e-04 2.696601e-04  5.653370e-06
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

6.3.2.3 Prediction plot for marginal effect of conductivity on median infection intensity with Trichodina

newdata1 <- newdata; newdata1[,"Con_av"] <- environment2$Con_av
BMA_Tri_Con_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3h = ggplot(environment2, aes(Con_av, BMA_Tri_Con_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_Tri_Con_av$fit-BMA_Tri_Con_av$se.bma.fit), ymax = (BMA_Tri_Con_av$fit+BMA_Tri_Con_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=Con_av, y=medin$Tri)) +
  labs(x=expression(paste("Conductivity [", mu, "S/cm]")), y=expression(paste(italic("Trichodina"), " sp. - infection intensity"))) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "h)")
figure3h

6.3.2.3 Prediction plot for marginal effect of COD on median infection intensity with Trichodina

newdata1 <- newdata; newdata1[,"COD_av"] <- environment2$COD_av
BMA_Tri_COD_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3d = ggplot(environment2, aes(COD_av, BMA_Tri_COD_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_Tri_COD_av$fit-BMA_Tri_COD_av$se.bma.fit), ymax = (BMA_Tri_COD_av$fit+BMA_Tri_COD_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=COD_av, y=medin$Tri)) +
  labs(x=expression(paste("Chemical Oxygen Demand [mg/L]")), y=expression(paste(italic("Trichodina"), " sp. - infection intensity"))) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "d)")
figure3d

6.4.1 Prevalence

bas.model <- bas.lm(prev$Tri ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)   model 1   model 2    model 3   model 4   model 5
## Intercept       1.00000000 1.0000000 1.0000000  1.0000000 1.0000000  1.000000
## avlength        0.04005815 0.0000000 0.0000000  0.0000000 0.0000000  0.000000
## avcondition     0.03389439 0.0000000 0.0000000  0.0000000 0.0000000  0.000000
## T_av            0.03654853 0.0000000 0.0000000  0.0000000 0.0000000  0.000000
## O2_sat_av       0.03749123 0.0000000 0.0000000  0.0000000 0.0000000  0.000000
## Con_av          0.20573471 0.0000000 1.0000000  0.0000000 1.0000000  0.000000
## COD_av          0.04629842 0.0000000 0.0000000  0.0000000 0.0000000  0.000000
## NH4._av         0.06690739 0.0000000 0.0000000  1.0000000 0.0000000  0.000000
## Nt_av           0.04604656 0.0000000 0.0000000  0.0000000 0.0000000  1.000000
## pool_riffle1    0.03411173 0.0000000 0.0000000  0.0000000 0.0000000  0.000000
## meander1        0.03657952 0.0000000 0.0000000  0.0000000 0.0000000  0.000000
## netcen          0.07467299 0.0000000 0.0000000  0.0000000 1.0000000  0.000000
## updist          0.07352673 0.0000000 0.0000000  0.0000000 0.0000000  0.000000
## BF                      NA 0.5616705 0.7274201  0.2732671 1.0000000  0.155205
## PostProbs               NA 0.6354000 0.0686000  0.0258000 0.0171000  0.014600
## R2                      NA 0.0000000 0.1264000  0.0750000 0.2249000  0.044000
## dim                     NA 1.0000000 2.0000000  2.0000000 3.0000000  2.000000
## logmarg                 NA 0.0000000 0.2585887 -0.7204656 0.5768398 -1.286169
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                       2.5%        97.5%          beta
## Intercept     7.818130e-01 8.857532e-01  8.340580e-01
## avlength      0.000000e+00 0.000000e+00 -1.453558e-04
## avcondition   0.000000e+00 0.000000e+00 -2.610917e-03
## T_av          0.000000e+00 0.000000e+00 -1.135202e-05
## O2_sat_av     0.000000e+00 0.000000e+00 -1.893466e-05
## Con_av       -1.787956e-07 4.202091e-04  6.570599e-05
## COD_av        0.000000e+00 0.000000e+00  9.451636e-05
## NH4._av      -3.379414e-05 1.957452e-02  1.713706e-03
## Nt_av         0.000000e+00 0.000000e+00  1.739937e-04
## pool_riffle1  0.000000e+00 0.000000e+00  2.104077e-04
## meander1      0.000000e+00 0.000000e+00 -6.652944e-04
## netcen       -3.742651e-06 9.178620e-09 -3.317353e-07
## updist       -1.636752e-06 5.369651e-10 -1.309710e-07
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
write.table(confint, 'GyroAA.txt', sep="\t")
pip <- summary(bas.model)
PIP[c(1:12),3] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

6.5 Variation in Glugea infection

bas.model <- bas.glm(pa$Glu ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av + 
                       NH4._av + Nt_av + SM_av + pool_riffle + meander + 
                       spavar$netcen + spavar$updist, data=environment2, betaprior=g.prior(100), family=binomial)
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##               P(B != 0 | Y)   model 1    model 2   model 3    model 4
## Intercept         1.0000000  1.000000  1.0000000  1.000000  1.0000000
## avlength          0.6721313  0.000000  1.0000000  0.000000  0.0000000
## avcondition       0.4472412  0.000000  1.0000000  1.000000  1.0000000
## T_av              0.2835083  0.000000  0.0000000  0.000000  0.0000000
## O2_sat_av         0.3836304  0.000000  0.0000000  0.000000  0.0000000
## Con_av            0.9992554  1.000000  1.0000000  1.000000  1.0000000
## COD_av            0.5410645  1.000000  1.0000000  0.000000  1.0000000
## NH4._av           0.8430908  1.000000  1.0000000  1.000000  1.0000000
## Nt_av             0.5500488  1.000000  0.0000000  1.000000  1.0000000
## SM_av             0.4324219  1.000000  0.0000000  1.000000  0.0000000
## pool_riffle1      0.4305176  0.000000  0.0000000  0.000000  0.0000000
## meander1          0.9996338  1.000000  1.0000000  1.000000  1.0000000
## spavar$netcen     0.4574951  0.000000  0.0000000  1.000000  0.0000000
## spavar$updist     0.8283325  1.000000  1.0000000  0.000000  1.0000000
## BF                       NA  1.000000  0.8189871  0.769804  0.7351301
## PostProbs                NA  0.053100  0.0522000  0.039900  0.0384000
## R2                       NA  1.000000  1.0000000  1.000000  1.0000000
## dim                      NA  8.000000  8.0000000  8.000000  8.0000000
## logmarg                  NA -5.808403 -6.0080902 -6.070023 -6.1161110
##                  model 5
## Intercept      1.0000000
## avlength       1.0000000
## avcondition    0.0000000
## T_av           0.0000000
## O2_sat_av      1.0000000
## Con_av         1.0000000
## COD_av         0.0000000
## NH4._av        0.0000000
## Nt_av          0.0000000
## SM_av          0.0000000
## pool_riffle1   1.0000000
## meander1       1.0000000
## spavar$netcen  1.0000000
## spavar$updist  1.0000000
## BF             0.5830791
## PostProbs      0.0315000
## R2             1.0000000
## dim            8.0000000
## logmarg       -6.3478357
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE
confint(coef.model)
##                        2.5%        97.5%          beta
## Intercept     -2.877566e+06 3.000550e+06 -6.056996e+03
## avlength      -4.420045e+04 4.258743e+04  8.562288e+01
## avcondition   -8.578579e+05 9.088818e+05 -1.218086e+01
## T_av          -3.140698e+04 2.899503e+04  1.760489e+01
## O2_sat_av     -5.868290e+03 5.525430e+03  8.380211e+00
## Con_av        -1.518537e+03 1.619110e+03  4.195005e+00
## COD_av        -1.323990e+04 1.518345e+04  9.690644e+00
## NH4._av       -2.252463e+05 2.166739e+05 -2.691592e+02
## Nt_av         -4.470842e+04 4.775553e+04  3.669994e+01
## SM_av         -2.069908e+03 2.430891e+03 -1.161863e+00
## pool_riffle1  -1.551102e+05 1.439102e+05  8.409678e+01
## meander1      -4.719167e+05 5.411526e+05 -1.301335e+03
## spavar$netcen -8.973198e+00 1.065294e+01 -8.483558e-03
## spavar$updist -5.016680e+00 5.104948e+00 -7.734368e-03
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),9] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

6.5 Variation in Contracaecum infection

bas.model <- bas.glm(pa$Con ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av + 
                       NH4._av + Nt_av + SM_av + pool_riffle + meander + 
                       spavar$netcen + spavar$updist, data=environment2, betaprior=g.prior(100), family=binomial)
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##               P(B != 0 | Y)   model 1     model 2     model 3     model 4
## Intercept        1.00000000   1.00000   1.0000000   1.0000000   1.0000000
## avlength         0.08128662   0.00000   1.0000000   0.0000000   0.0000000
## avcondition      0.02805176   0.00000   0.0000000   1.0000000   0.0000000
## T_av             0.01058350   0.00000   0.0000000   0.0000000   0.0000000
## O2_sat_av        0.02976074   0.00000   0.0000000   0.0000000   1.0000000
## Con_av           0.01105957   0.00000   0.0000000   0.0000000   0.0000000
## COD_av           0.01395264   0.00000   0.0000000   0.0000000   0.0000000
## NH4._av          0.01990967   0.00000   0.0000000   0.0000000   0.0000000
## Nt_av            0.01219482   0.00000   0.0000000   0.0000000   0.0000000
## SM_av            0.01367188   0.00000   0.0000000   0.0000000   0.0000000
## pool_riffle1     0.01322021   0.00000   0.0000000   0.0000000   0.0000000
## meander1         0.01164551   0.00000   0.0000000   0.0000000   0.0000000
## spavar$netcen    0.01945801   0.00000   0.0000000   0.0000000   0.0000000
## spavar$updist    0.01293945   0.00000   0.0000000   0.0000000   0.0000000
## BF                       NA   1.00000   0.9507459   0.2581488   0.2583995
## PostProbs                NA   0.77910   0.0562000   0.0172000   0.0172000
## R2                       NA   0.00000   0.0864000   0.0365000   0.0366000
## dim                      NA   1.00000   2.0000000   2.0000000   2.0000000
## logmarg                  NA -25.82594 -25.8764468 -27.1801577 -27.1791867
##                 model 5
## Intercept       1.00000
## avlength        0.00000
## avcondition     0.00000
## T_av            0.00000
## O2_sat_av       0.00000
## Con_av          0.00000
## COD_av          0.00000
## NH4._av         1.00000
## Nt_av           0.00000
## SM_av           0.00000
## pool_riffle1    0.00000
## meander1        0.00000
## spavar$netcen   0.00000
## spavar$updist   0.00000
## BF              0.19791
## PostProbs       0.01350
## R2              0.02640
## dim             2.00000
## logmarg       -27.44588
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE
confint(coef.model)
##                    2.5%     97.5%          beta
## Intercept     -9.880859 2.3742736 -5.430024e-01
## avlength       0.000000 0.1497804  1.509881e-02
## avcondition    0.000000 0.0000000 -1.814762e-01
## T_av           0.000000 0.0000000  7.945717e-04
## O2_sat_av      0.000000 0.0000000  1.286763e-03
## Con_av         0.000000 0.0000000  5.456325e-06
## COD_av         0.000000 0.0000000 -3.836691e-04
## NH4._av        0.000000 0.0000000 -6.098874e-03
## Nt_av          0.000000 0.0000000 -8.846247e-04
## SM_av          0.000000 0.0000000  1.337787e-04
## pool_riffle1   0.000000 0.0000000 -5.842378e-03
## meander1       0.000000 0.0000000  3.379776e-03
## spavar$netcen  0.000000 0.0000000 -8.679200e-07
## spavar$updist  0.000000 0.0000000 -1.196324e-07
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),10] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

Variation in Anguillicoloides infection

bas.model <- bas.glm(pa$Ang ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av + 
                       NH4._av + Nt_av + SM_av + pool_riffle + meander + 
                       spavar$netcen + spavar$updist, data=environment2, betaprior=g.prior(100), family=binomial)
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##               P(B != 0 | Y)   model 1     model 2     model 3     model 4
## Intercept         1.0000000  1.000000   1.0000000   1.0000000   1.0000000
## avlength          0.9849854  1.000000   1.0000000   1.0000000   1.0000000
## avcondition       0.3566895  0.000000   0.0000000   0.0000000   0.0000000
## T_av              0.9871948  1.000000   1.0000000   1.0000000   1.0000000
## O2_sat_av         0.8635742  1.000000   1.0000000   1.0000000   1.0000000
## Con_av            0.5832520  0.000000   1.0000000   1.0000000   0.0000000
## COD_av            0.3989502  0.000000   0.0000000   1.0000000   0.0000000
## NH4._av           0.9910645  1.000000   1.0000000   1.0000000   1.0000000
## Nt_av             0.5583984  0.000000   1.0000000   0.0000000   0.0000000
## SM_av             0.9698608  1.000000   1.0000000   1.0000000   1.0000000
## pool_riffle1      0.4013184  0.000000   0.0000000   0.0000000   1.0000000
## meander1          0.9865601  1.000000   1.0000000   1.0000000   1.0000000
## spavar$netcen     0.8922974  1.000000   1.0000000   1.0000000   1.0000000
## spavar$updist     0.5290039  1.000000   0.0000000   0.0000000   1.0000000
## BF                       NA  1.000000   0.4728289   0.3794633   0.2858158
## PostProbs                NA  0.080800   0.0735000   0.0537000   0.0429000
## R2                       NA  1.000000   1.0000000   1.0000000   1.0000000
## dim                      NA  9.000000  10.0000000  10.0000000  10.0000000
## logmarg                  NA -9.802499 -10.5515211 -10.7714968 -11.0549070
##                   model 5
## Intercept       1.0000000
## avlength        1.0000000
## avcondition     1.0000000
## T_av            1.0000000
## O2_sat_av       0.0000000
## Con_av          1.0000000
## COD_av          0.0000000
## NH4._av         1.0000000
## Nt_av           1.0000000
## SM_av           1.0000000
## pool_riffle1    0.0000000
## meander1        1.0000000
## spavar$netcen   1.0000000
## spavar$updist   0.0000000
## BF              0.2037316
## PostProbs       0.0324000
## R2              1.0000000
## dim            10.0000000
## logmarg       -11.3934513
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE
confint(coef.model)
##                        2.5%        97.5%          beta
## Intercept     -2.427843e+06 2.407495e+06 -6.443834e+03
## avlength      -4.661497e+04 4.211047e+04  1.321190e+02
## avcondition   -8.356648e+05 8.401762e+05 -3.035671e+02
## T_av          -1.531048e+05 1.561560e+05  5.865874e+02
## O2_sat_av     -1.899200e+04 1.821646e+04 -3.714070e+01
## Con_av        -6.375893e+02 7.129852e+02  3.717928e-01
## COD_av        -7.961092e+03 7.947834e+03  1.910391e+00
## NH4._av       -3.591783e+05 3.383641e+05 -1.053055e+03
## Nt_av         -4.799272e+04 3.935074e+04 -4.605084e+01
## SM_av         -5.743284e+03 5.697560e+03  1.280831e+01
## pool_riffle1  -1.496669e+05 1.333208e+05 -1.219107e+02
## meander1      -4.050513e+05 3.841885e+05 -1.391517e+03
## spavar$netcen -1.550874e+01 1.393704e+01 -2.782679e-02
## spavar$updist -3.150044e+00 2.889053e+00  1.171764e-03
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),11] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

Variation in Individual Parasitation Index (all parasites)

bas.model <- bas.lm(avPI ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)    model 1    model 2    model 3  model 4  model 5
## Intercept        1.0000000 1.00000000 1.00000000 1.00000000 1.000000 1.000000
## avlength         0.2261960 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## avcondition      0.1810890 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## T_av             0.4574073 0.00000000 1.00000000 0.00000000 1.000000 0.000000
## O2_sat_av        0.1952767 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## Con_av           0.4598919 0.00000000 0.00000000 0.00000000 0.000000 1.000000
## COD_av           0.5684747 0.00000000 0.00000000 0.00000000 0.000000 1.000000
## NH4._av          0.2748144 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## Nt_av            0.4582802 0.00000000 0.00000000 1.00000000 1.000000 0.000000
## pool_riffle1     0.3054914 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## meander1         0.1985617 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## netcen           0.1910872 0.00000000 0.00000000 0.00000000 0.000000 0.000000
## updist           0.4153657 0.00000000 0.00000000 0.00000000 0.000000 1.000000
## BF                      NA 0.01969453 0.09251289 0.08737546 0.355023 1.000000
## PostProbs               NA 0.10850000 0.04250000 0.04010000 0.029600 0.025000
## R2                      NA 0.00000000 0.18980000 0.18710000 0.325400 0.424800
## dim                     NA 1.00000000 2.00000000 2.00000000 3.000000 4.000000
## logmarg                 NA 0.00000000 1.54700725 1.48987379 2.891842 3.927415
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                       2.5%        97.5%          beta
## Intercept     1.116139e+00 1.708031e+00  1.422090e+00
## avlength     -8.564190e-02 1.467677e-02 -7.992245e-03
## avcondition  -2.839253e+00 1.768409e+00 -9.325811e-02
## T_av         -2.523565e-04 3.639709e-01  8.644408e-02
## O2_sat_av    -2.287262e-02 1.596544e-02 -7.476618e-04
## Con_av       -6.334651e-07 3.201875e-03  8.405554e-04
## COD_av        0.000000e+00 6.203673e-02  2.109188e-02
## NH4._av      -4.940566e-01 1.660406e-01 -4.965441e-02
## Nt_av         0.000000e+00 2.277704e-01  5.507750e-02
## pool_riffle1 -4.507692e-04 9.323231e-01  1.389140e-01
## meander1     -5.556519e-01 2.496543e-01 -3.866036e-02
## netcen       -2.956977e-05 2.000712e-05 -1.270941e-06
## updist       -2.688456e-05 9.716094e-09 -5.952452e-06
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),12] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

Prediction plot for marginal effect of COD on Individual Parasitation Index

newdata1 <- newdata; newdata1[,"COD_av"] <- environment2$COD_av
BMA_IPI_COD_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3b = ggplot(environment2, aes(COD_av, BMA_IPI_COD_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_IPI_COD_av$fit-BMA_IPI_COD_av$se.bma.fit), ymax = (BMA_IPI_COD_av$fit+BMA_IPI_COD_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=COD_av, y=avPI)) +
  labs(x=expression("Chemical Oxygen Demand [mg/L]"), y=expression("I"[PI]*" - all parasites")) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "b)")
figure3b

Variation in Individual Parasitation Index (only ectoparasites)

bas.model <- bas.lm(avPI_ecto ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)   model 1      model 2   model 3  model 4  model 5
## Intercept        1.0000000 1.0000000  1.000000000 1.0000000 1.000000 1.000000
## avlength         0.4669536 0.0000000  1.000000000 0.0000000 0.000000 1.000000
## avcondition      0.3758074 0.0000000  1.000000000 0.0000000 0.000000 0.000000
## T_av             0.2753561 0.0000000  1.000000000 0.0000000 0.000000 0.000000
## O2_sat_av        0.2919614 0.0000000  1.000000000 0.0000000 0.000000 0.000000
## Con_av           0.6617062 0.0000000  1.000000000 1.0000000 0.000000 0.000000
## COD_av           0.8437205 0.0000000  1.000000000 1.0000000 1.000000 1.000000
## NH4._av          0.4163705 0.0000000  1.000000000 0.0000000 0.000000 0.000000
## Nt_av            0.8302866 1.0000000  1.000000000 0.0000000 1.000000 1.000000
## pool_riffle1     0.4025108 0.0000000  1.000000000 0.0000000 0.000000 0.000000
## meander1         0.5031625 0.0000000  1.000000000 0.0000000 0.000000 0.000000
## netcen           0.3010717 0.0000000  1.000000000 0.0000000 0.000000 0.000000
## updist           0.3285744 0.0000000  1.000000000 0.0000000 0.000000 0.000000
## BF                      NA 0.1318791  0.008313608 0.4282678 0.319205 1.000000
## PostProbs               NA 0.0510000  0.038600000 0.0301000 0.022500 0.021100
## R2                      NA 0.3007000  0.682800000 0.4141000 0.403700 0.496600
## dim                     NA 2.0000000 13.000000000 3.0000000 3.000000 4.000000
## logmarg                 NA 4.0635523  1.299560343 5.2414154 4.947500 6.089422
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                       2.5%        97.5%          beta
## Intercept     6.899850e-01 1.159346e+00  9.199577e-01
## avlength     -1.072538e-01 6.599823e-03 -2.357999e-02
## avcondition  -4.610896e+00 3.985984e-01 -7.055549e-01
## T_av         -9.506040e-02 1.406152e-01  5.818763e-03
## O2_sat_av    -1.638007e-02 2.228484e-02  1.122707e-03
## Con_av        0.000000e+00 3.193247e-03  1.203349e-03
## COD_av        0.000000e+00 6.428436e-02  3.405818e-02
## NH4._av      -4.990487e-01 7.994328e-02 -8.034086e-02
## Nt_av        -3.936480e-05 2.546442e-01  1.257830e-01
## pool_riffle1 -7.124275e-02 8.162781e-01  1.490513e-01
## meander1     -9.338322e-01 2.458840e-02 -2.358742e-01
## netcen       -3.428514e-05 2.027454e-05 -2.447510e-06
## updist       -1.553765e-05 4.040408e-06 -1.836166e-06
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))

## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),13] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

Prediction plot for marginal effect of COD on Individual Parasitation Index (only ectoparasites)

newdata1 <- newdata; newdata1[,"COD_av"] <- environment2$COD_av
BMA_IPIecto_COD_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3c = ggplot(environment2, aes(COD_av, BMA_IPIecto_COD_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_IPIecto_COD_av$fit-BMA_IPIecto_COD_av$se.bma.fit), ymax = (BMA_IPIecto_COD_av$fit+BMA_IPIecto_COD_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=COD_av, y=avPI_ecto)) +
  labs(x=expression("Chemical Oxygen Demand [mg/L]"), y=expression("I"[PI]*" - ectoparasites")) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "c)")
figure3c

Prediction plot for marginal effect of total nitrogen on Individual Parasitation Index (only ectoparasites)

newdata1 <- newdata; newdata1[,"Nt_av"] <- environment2$Nt_av
BMA_IPIecto_Nt_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3f = ggplot(environment2, aes(Nt_av, BMA_IPIecto_Nt_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_IPIecto_Nt_av$fit-BMA_IPIecto_Nt_av$se.bma.fit), ymax = (BMA_IPIecto_Nt_av$fit+BMA_IPIecto_Nt_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=Nt_av, y=avPI_ecto)) +
  labs(x=expression("Nt [mg/L]"), y=expression("I"[PI]*" - ectoparasites")) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "f)")
figure3f

Prediction plot for marginal effect of conductivity on Individual Parasitation Index (only ectoparasites)

newdata1 <- newdata; newdata1[,"Con_av"] <- environment2$Con_av
BMA_IPIecto_Con_av <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3g = ggplot(environment2, aes(Con_av, BMA_IPIecto_Con_av$fit)) +
  theme_bw() +
  geom_line(color="red", size=1) +
  geom_ribbon(aes(ymin = (BMA_IPIecto_Con_av$fit-BMA_IPIecto_Con_av$se.bma.fit), ymax = (BMA_IPIecto_Con_av$fit+BMA_IPIecto_Con_av$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=Con_av, y=avPI)) +
  labs(x=expression(paste("Conductivity [", mu, "S/cm]")), y=expression("I"[PI]*" - all parasites")) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "g)")
figure3g

Prediction plot for marginal effect of meanders on Individual Parasitation Index (only ectoparasites)

newdata1 <- newdata; newdata1[,"meander"] <- environment2$meander
BMA_IPIecto_meander <- predict(bas.model, newdata = newdata1, estimator = "BMA", se.fit=TRUE)
figure3l = ggplot(environment2, aes(meander, BMA_IPIecto_meander$fit)) +
  theme_bw() +
  #geom_line(color="red", size=1) +
  #geom_ribbon(aes(ymin = (BMA$fit-BMA$se.bma.fit), ymax = (BMA$fit+BMA$se.bma.fit)), alpha = .1) +
  geom_point(data = environment2, aes(x=meander, y=avPI)) +
  geom_boxplot(aes(lower = (BMA_IPIecto_meander$fit-BMA_IPIecto_meander$se.bma.fit), middle = BMA_IPIecto_meander$fit, upper = (BMA_IPIecto_meander$fit+BMA_IPIecto_meander$se.bma.fit))) +
  labs(x=expression("Meander"), y=expression("I"[PI]*" - ectoparasites")) +
  theme(axis.title.x = element_text(size=12),
        axis.title.y = element_text(size=12)) +  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  labs(subtitle = "l)")
figure3l

Variation in Individual Parasitation Index (only endoparasites)

bas.model <- bas.lm(avPI_endo ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2, prior="JZS")
yhat = fitted(bas.model, estimator = "BMA") #these are the fitted values under BMA
r = bas.model$Y - yhat #these are the model residuals
plot(bas.model)

summary(bas.model)
##              P(B != 0 | Y)   model 1  model 2     model 3    model 4   model 5
## Intercept       1.00000000 1.0000000 1.000000  1.00000000  1.0000000 1.0000000
## avlength        0.06442140 0.0000000 0.000000  0.00000000  0.0000000 0.0000000
## avcondition     0.12246410 0.0000000 0.000000  1.00000000  0.0000000 0.0000000
## T_av            0.30397004 0.0000000 1.000000  0.00000000  0.0000000 1.0000000
## O2_sat_av       0.05178919 0.0000000 0.000000  0.00000000  0.0000000 0.0000000
## Con_av          0.05158494 0.0000000 0.000000  0.00000000  0.0000000 0.0000000
## COD_av          0.04593088 0.0000000 0.000000  0.00000000  0.0000000 0.0000000
## NH4._av         0.08528364 0.0000000 0.000000  0.00000000  0.0000000 1.0000000
## Nt_av           0.06477759 0.0000000 0.000000  0.00000000  0.0000000 0.0000000
## pool_riffle1    0.04504579 0.0000000 0.000000  0.00000000  0.0000000 0.0000000
## meander1        0.05899866 0.0000000 0.000000  0.00000000  0.0000000 0.0000000
## netcen          0.05414430 0.0000000 0.000000  0.00000000  0.0000000 0.0000000
## updist          0.10936917 0.0000000 0.000000  0.00000000  1.0000000 0.0000000
## BF                      NA 0.3484238 1.000000  0.33059983  0.2877395 0.6134470
## PostProbs               NA 0.4921000 0.117700  0.03890000  0.0339000 0.0131000
## R2                      NA 0.0000000 0.166100  0.11040000  0.1032000 0.2244000
## dim                     NA 1.0000000 2.000000  2.00000000  2.0000000 3.0000000
## logmarg                 NA 0.0000000 1.054336 -0.05251095 -0.1913639 0.5656742
image(bas.model, rotate=F)

coef.model <- coef(bas.model)
abs(coef.model$postmean)-2*coef.model$postsd > 0
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE
confint(coef.model)
##                       2.5%        97.5%          beta
## Intercept     2.974020e-01 6.917548e-01  4.947994e-01
## avlength     -3.456166e-04 1.669276e-02  1.398963e-03
## avcondition  -3.455674e-03 2.596667e+00  2.446911e-01
## T_av          0.000000e+00 2.137081e-01  4.520999e-02
## O2_sat_av    -3.702737e-04 3.813693e-05 -2.446396e-04
## Con_av       -6.714701e-05 4.289287e-06 -1.385311e-05
## COD_av        0.000000e+00 0.000000e+00 -1.072193e-04
## NH4._av      -1.119913e-01 1.101911e-04 -9.322620e-03
## Nt_av        -2.836973e-02 1.844054e-03 -2.149910e-03
## pool_riffle1  0.000000e+00 0.000000e+00  2.267863e-03
## meander1     -1.340521e-02 7.816289e-02  1.005422e-02
## netcen       -1.286967e-06 8.497622e-08 -3.607730e-07
## updist       -9.420533e-06 0.000000e+00 -8.091451e-07
## attr(,"Probability")
## [1] 0.95
## attr(,"class")
## [1] "confint.bas"
plot(confint(coef.model, parm = 2:11))
## Warning in arrows(x[not.deg], ci[not.deg, 1], x[not.deg], ci[not.deg, 2], :
## zero-length arrow is of indeterminate angle and so skipped

## NULL
confint <- confint(coef.model, parm = 2:11)
pip <- summary(bas.model)
PIP[c(1:12),14] <- pip[2:13,1]*sign(coef.model$postmean[2:13])

Figure 3

library(gridExtra)
pdf("Figure3.pdf", width = 10.8, height = 14.4)
grid.arrange(figure3a, figure3b, figure3c, figure3d, figure3e, figure3f, figure3g, figure3h, figure3i, figure3j, figure3k, figure3l)
dev.off()
## png 
##   2

7. BORAL analysis

Model-based analysis of multivariate abundance data using Bayesian Markov chain Monte Carlo methods for parameter estimation

7.1 BORAL analysis for average abundances of parasites

data$Site <- as.factor(data$site)
levels(data$site) <- levels(as.factor(environment2$Site))
data_m <- merge(data, environment2, by = "Site")
data_all <- na.omit(data_m) 
names(data_all)
##  [1] "Site"                "site"                "fish"               
##  [4] "parspeciesrichness"  "div_shannon"         "div_simpson"        
##  [7] "temp"                "pH"                  "conductivity"       
## [10] "nitrogen"            "phosphorus"          "oxygen"             
## [13] "netcen.x"            "updist.x"            "updist2"            
## [16] "updist3"             "fishspeciesrichness" "weight"             
## [19] "weigh..g."           "length"              "SMI"                
## [22] "Sex"                 "Gyr"                 "Tri"                
## [25] "Glu"                 "ecto_screener"       "Con"                
## [28] "CysL"                "Pro"                 "Aca"                
## [31] "Cam"                 "Ang"                 "CysI"               
## [34] "endo_screener"       "PI"                  "PI_ecto"            
## [37] "PI_endo"             "T_av"                "O2_sat_av"          
## [40] "Con_av"              "COD_av"              "NH4._av"            
## [43] "Nt_av"               "SM_av"               "pool_riffle"        
## [46] "meander"             "updist.y"            "netcen.y"
avcondition <- aggregate(data$SMI, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]
avlength <- aggregate(data$length, by = list(data[,1]), function(x){mean(x, na.rm =T)})[,2]

y <- round(cbind(avab$Gyr, avab$Tri, avab$Glu, avab$Con, avab$Ang))
X <- cbind(avcondition,
           avlength,
           environment2$T_av,
           environment2$O2_sat_av,
           environment2$Con_av,
           environment2$COD_av, 
           environment2$NH4._av, 
           environment2$Nt_av, 
           environment2$netcen, 
           environment2$updist, 
           as.numeric(environment2$pool_riffle), 
           as.numeric(environment2$meander))
colnames(X) <- c("avcondition", "avlength", "T", "O2", "Con", "COD", "NH4", "Nt", "netcen", "updist", "pool_riffle", "meander")

example_mcmc_control <- list(n.burnin = 1000, n.iteration = 10000, n.thin = 1)
testpath <- file.path(tempdir(), "jagsboralmodel.txt")
paramod <- boral(y, X = X,
                      family = "negative.binomial",
                      mcmc.control = example_mcmc_control,
                      model.name = testpath,
                      lv.control = list(num.lv = 2, type = "independent"),
                      save.model = TRUE)
## module glm loaded
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 185
##    Unobserved stochastic nodes: 338
##    Total graph size: 2173
## 
## Initializing model
plot(paramod)
## NULL

coefsplot(covname = "avcondition", object = paramod) #Condition

coefsplot(covname = "avlength", object = paramod) #Length

coefsplot(covname = "T", object = paramod) #Temperature

coefsplot(covname = "O2", object = paramod) #Oxygen

coefsplot(covname = "Con", object = paramod) #Conductivity

coefsplot(covname = "COD", object = paramod) #COD

coefsplot(covname = "NH4", object = paramod) #NH4

coefsplot(covname = "Nt", object = paramod) #Nt

coefsplot(covname = "netcen", object = paramod) #netcen

coefsplot(covname = "updist", object = paramod) #updist

coefsplot(covname = "pool_riffle", object = paramod) #poolriffle

coefsplot(covname = "meander", object = paramod) #meander

envcors <- get.enviro.cor(paramod)
rescors <- get.residual.cor(paramod)
corrplot(envcors$sig.cor, type = "lower", diag = FALSE, title = "Correlations due to covariates", mar = c(3,0.5,2,1), tl.srt = 45)

corrplot(rescors$sig.cor, type = "lower", diag = FALSE, title = "Residual correlations", mar = c(3,0.5,2,1), tl.srt = 45)

7.1 BORAL analysis for median infection intensities of parasites

y <- round(cbind(medin$Gyr, medin$Tri, medin$Glu, medin$Con, medin$Ang))
paramod <- boral(y, X = X,
                      family = "negative.binomial",
                      mcmc.control = example_mcmc_control,
                      model.name = testpath,
                      lv.control = list(num.lv = 2, type = "independent"),
                      save.model = TRUE)
## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 185
##    Unobserved stochastic nodes: 338
##    Total graph size: 2173
## 
## Initializing model
plot(paramod)
## NULL

coefsplot(covname = "avcondition", object = paramod) #Condition

coefsplot(covname = "avlength", object = paramod) #Length

coefsplot(covname = "T", object = paramod) #Temperature

coefsplot(covname = "O2", object = paramod) #Oxygen

coefsplot(covname = "Con", object = paramod) #Conductivity

coefsplot(covname = "COD", object = paramod) #COD

coefsplot(covname = "NH4", object = paramod) #NH4

coefsplot(covname = "Nt", object = paramod) #Nt

coefsplot(covname = "netcen", object = paramod) #netcen

coefsplot(covname = "updist", object = paramod) #updist

coefsplot(covname = "pool_riffle", object = paramod) #poolriffle

coefsplot(covname = "meander", object = paramod) #meander

envcors <- get.enviro.cor(paramod)
rescors <- get.residual.cor(paramod)
corrplot(envcors$sig.cor, type = "lower", diag = FALSE, title = "Correlations due to covariates", mar = c(3,0.5,2,1), tl.srt = 45)

corrplot(rescors$sig.cor, type = "lower", diag = FALSE, title = "Residual correlations", mar = c(3,0.5,2,1), tl.srt = 45)

8. Multivariate analysis

8.1 Component communities

# Component communities: Bray-Curtis dissimilarities based on Hellinger transformed average abundance data
spe.hel_bray <- vegdist(decostand(avab[,-1], na.rm=T, method="hellinger"), method="bray", na.rm=T)

# Check whether Euclidean and Bray-Curtis distances are comparable
spe.hel_euc <- vegdist(decostand(avab[,-1], na.rm=T, method="hellinger"), method="euc", na.rm=T)
plot(spe.hel_bray, spe.hel_euc)

mantel(spe.hel_bray, spe.hel_euc)
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = spe.hel_bray, ydis = spe.hel_euc) 
## 
## Mantel statistic r: 0.9648 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.105 0.138 0.165 0.196 
## Permutation: free
## Number of permutations: 999
model_adonis = adonis(spe.hel_bray ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2)
## 'adonis' will be deprecated: use 'adonis2' instead
model_adonis$aov.tab
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##             Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## avlength     1    0.1863 0.18629  1.8398 0.04303  0.117  
## avcondition  1    0.1305 0.13055  1.2893 0.03016  0.267  
## T_av         1    0.3204 0.32039  3.1643 0.07401  0.016 *
## O2_sat_av    1    0.1368 0.13678  1.3509 0.03160  0.241  
## Con_av       1    0.1326 0.13262  1.3098 0.03064  0.266  
## COD_av       1    0.0657 0.06568  0.6486 0.01517  0.668  
## NH4._av      1    0.2040 0.20399  2.0146 0.04712  0.096 .
## Nt_av        1    0.1451 0.14509  1.4329 0.03352  0.227  
## pool_riffle  1    0.0853 0.08529  0.8424 0.01970  0.517  
## meander      1    0.2029 0.20286  2.0035 0.04686  0.098 .
## netcen       1    0.2173 0.21728  2.1459 0.05019  0.088 .
## updist       1    0.0719 0.07193  0.7104 0.01662  0.621  
## Residuals   24    2.4301 0.10125         0.56137         
## Total       36    4.3288                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# environmental variables
env_select <- environment2[,c("T_av", "O2_sat_av", "Con_av", "COD_av", "NH4._av", "Nt_av", "pool_riffle", "meander", "netcen", "updist")]
env_select$pool_riffle <- as.numeric(env_select$pool_riffle)
env_select$meander <- as.numeric(env_select$meander)

pca <- prcomp(env_select, scale.=T)
summary(pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     1.7124 1.5545 1.1221 1.0140 0.88807 0.79463 0.56647
## Proportion of Variance 0.2933 0.2416 0.1259 0.1028 0.07887 0.06314 0.03209
## Cumulative Proportion  0.2933 0.5349 0.6608 0.7636 0.84248 0.90563 0.93771
##                            PC8     PC9    PC10
## Standard deviation     0.50483 0.46939 0.38429
## Proportion of Variance 0.02549 0.02203 0.01477
## Cumulative Proportion  0.96320 0.98523 1.00000
plot(pca)

biplot(pca)

# Assess the effect of environmental variables on parasite component community dissimilarities using distance based RDA
spe.rda <- dbrda(spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
summary(spe.rda)
## 
## Call:
## dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av +      NH4._av + Nt_av + pool_riffle + meander, data = env_select) 
## 
## Partitioning of squared Bray distance:
##               Inertia Proportion
## Total           4.329     1.0000
## Constrained     1.166     0.2695
## Unconstrained   3.162     0.7305
## 
## Eigenvalues, and their contribution to the squared Bray distance 
## 
## Importance of components:
##                       dbRDA1 dbRDA2  dbRDA3  dbRDA4   dbRDA5    dbRDA6
## Eigenvalue            0.5740 0.3372 0.17566 0.09415 0.042525 0.0020670
## Proportion Explained  0.1326 0.0779 0.04058 0.02175 0.009824 0.0004775
## Cumulative Proportion     NA     NA      NA      NA       NA        NA
##                         idbRDA1  idbRDA2   MDS1   MDS2   MDS3    MDS4    MDS5
## Eigenvalue            -0.014797 -0.04441 1.3249 0.8027 0.4680 0.31323 0.29121
## Proportion Explained   0.003418  0.01026 0.3061 0.1854 0.1081 0.07236 0.06727
## Cumulative Proportion        NA       NA     NA     NA     NA      NA      NA
##                          MDS6    MDS7    MDS8    MDS9   MDS10   MDS11   MDS12
## Eigenvalue            0.13644 0.11438 0.08987 0.07065 0.06425 0.02766 0.01576
## Proportion Explained  0.03152 0.02642 0.02076 0.01632 0.01484 0.00639 0.00364
## Cumulative Proportion      NA      NA      NA      NA      NA      NA      NA
##                          MDS13     MDS14      iMDS1     iMDS2    iMDS3
## Eigenvalue            0.011872 2.725e-04 -0.0020314 -0.007356 -0.01078
## Proportion Explained  0.002743 6.295e-05  0.0004693  0.001699  0.00249
## Cumulative Proportion       NA        NA         NA        NA       NA
##                           iMDS4    iMDS5     iMDS6     iMDS7    iMDS8     iMDS9
## Eigenvalue            -0.014744 -0.01935 -0.022495 -0.029849 -0.03342 -0.041489
## Proportion Explained   0.003406  0.00447  0.005197  0.006895  0.00772  0.009584
## Cumulative Proportion        NA       NA        NA        NA       NA        NA
##                         iMDS10   iMDS11   iMDS12   iMDS13   iMDS14
## Eigenvalue            -0.05984 -0.06439 -0.06995 -0.09116 -0.10203
## Proportion Explained   0.01382  0.01487  0.01616  0.02106  0.02357
## Cumulative Proportion       NA       NA       NA       NA       NA
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                       dbRDA1 dbRDA2 dbRDA3  dbRDA4  dbRDA5   dbRDA6  idbRDA1
## Eigenvalue            0.5740 0.3372 0.1757 0.09415 0.04252 0.002067 -0.01480
## Proportion Explained  0.4921 0.2891 0.1506 0.08072 0.03646 0.001772  0.01269
## Cumulative Proportion     NA     NA     NA      NA      NA       NA       NA
##                        idbRDA2
## Eigenvalue            -0.04441
## Proportion Explained   0.03807
## Cumulative Proportion       NA
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  3.533199 
## 
## 
## Site scores (weighted sums of species scores)
## 
##      dbRDA1   dbRDA2   dbRDA3    dbRDA4   dbRDA5   dbRDA6
## 1  -0.07530 -0.62580 -0.63602 -2.643196 -0.34427 -23.8946
## 2   0.87242  0.77832 -0.32068  2.074749 -0.25800  41.3754
## 3  -0.81094  1.29200 -0.57735 -0.708420 -2.35800  -5.0035
## 4  -0.97450  0.91051  1.79849 -2.192492 -1.99248  -4.7134
## 5  -0.25889  0.74559 -1.30995 -1.418227  0.33185 -21.9720
## 6   0.41640  1.01888 -0.34551 -2.318324  1.61924 -24.8949
## 7  -0.98286 -0.67495  0.62017  1.174218  1.55155  -0.8974
## 8   0.93807  0.62635 -3.10878 -0.765073  2.41981  -3.0764
## 9  -0.62816 -0.77245 -1.30540  1.806170  1.30543 -41.2679
## 10 -1.03727 -0.50089  0.62572  0.369727  1.00378   4.7819
## 11  1.10743 -0.99448 -0.66515 -0.071830  0.83402 -14.8528
## 12  0.61624  0.76828  0.87110  0.357836  0.12125  38.2552
## 13 -0.96088 -0.72453  0.47977  1.569869  1.43399  -6.8692
## 14  1.41128  2.03453 -1.67256  0.776255  0.27648  29.9435
## 15 -1.00448 -0.59302  0.65010  0.833950  1.44643   3.4902
## 16  1.26382  0.81623  3.19230  1.380128  1.12523 -10.4959
## 17 -0.92726 -0.73637  0.33774  1.892276  1.07787  -8.2400
## 18 -1.04523 -1.01280 -0.59052  2.897871  0.35870 -14.2865
## 19  0.89987 -1.54255 -0.54435  1.738188 -2.61939  32.8867
## 20 -0.14608 -0.34079  0.78630  1.316880  1.70249  -9.1246
## 21 -1.30124 -1.46318 -1.27370  3.739798  0.91928 -16.7232
## 22 -0.07485  1.73636  0.30728 -1.760677 -0.46041   6.1504
## 23 -0.18008 -0.64105  0.55428 -1.310793 -0.82437  -3.6579
## 24  0.53054  0.08088  1.31623 -3.529128 -1.96164  41.6937
## 25 -0.50901 -0.82222 -1.26480  2.167458 -6.71178 -22.4535
## 26  0.54874  0.92050  1.12409  0.009258  0.07714  34.0773
## 27  1.25804 -0.69768  2.37255 -1.697654 -3.24936 -16.2349
## 28 -0.86554  0.32492 -0.68347  0.704689 -0.16329  13.0277
## 29  1.82276  0.98120  0.33891 -1.852448  4.33581  -8.6587
## 30  0.72974 -0.59907 -1.00943 -0.201449  1.77124  -4.5186
## 31 -1.02904  0.19515  1.05102 -1.468951 -1.36610  23.1359
## 32 -0.13877 -0.06050 -0.19569 -1.899083 -1.10403 -31.4856
## 33  1.01209 -0.59332  2.58400 -0.318337 -1.43566  39.0635
## 34 -0.39319 -1.79743  0.20497 -0.362972 -0.59693 -19.8004
## 35 -0.89506  1.46497 -0.05622 -1.876722 -1.86120   3.5844
## 36  0.87386  1.17335 -3.50567  1.196922  1.51036  -9.4970
## 37 -0.06267 -0.67495 -0.14975  0.389535  2.08495  11.1529
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##      dbRDA1   dbRDA2   dbRDA3   dbRDA4    dbRDA5   dbRDA6
## 1  -0.16605 -0.69736 -0.45052  0.39949 -1.140870  0.26766
## 2  -0.07683 -0.33692  0.61638  0.50606 -0.381728  1.66676
## 3  -0.10071  0.76299  0.51759  0.02910  0.003366 -0.26259
## 4  -1.55548  0.90013  0.66933 -0.51241 -0.415318 -0.13093
## 5  -0.41914  0.27994 -0.66104 -1.10909  0.384267  0.50060
## 6   0.63090  0.30135 -0.07724  0.37780  0.857007 -0.33310
## 7  -0.15252 -0.25647 -0.20458 -0.38652  0.792931  0.48033
## 8  -0.14874  0.12938 -0.39866 -0.78891  0.083814  0.27677
## 9  -0.14359  0.11628  0.95439 -0.04629  0.609992 -0.64671
## 10 -0.74154 -0.65088  0.87072 -1.17151 -0.875610 -1.73263
## 11  0.77571 -0.48314  0.02321 -0.76003  0.259283 -0.01494
## 12 -0.44387 -0.28442  0.54938 -0.01813  1.062736 -0.71676
## 13 -0.25608 -0.61494 -0.47847  0.74269  0.180828 -0.14346
## 14 -0.23636  0.67472 -0.87185  0.72759  0.046993 -0.61306
## 15 -0.75508  0.19781  0.41035  0.89449 -0.744271 -0.25040
## 16  0.97946  1.50039  0.74911  0.12808 -0.657115 -0.30363
## 17 -1.19679 -0.18091 -0.90174  0.16093 -0.229189  0.70454
## 18  0.23197  0.39183 -0.42932 -0.82808 -0.242760  0.64991
## 19  0.80868 -0.55920 -0.22435  0.32426 -0.445610 -0.51746
## 20 -0.19206  0.12003  0.24728  0.99541  0.747110  0.40576
## 21 -0.38477  0.09371 -0.03632  1.07969  1.383306 -0.25819
## 22  0.02352  0.34014  0.07462  0.03662 -0.352806  0.35626
## 23  0.18716 -0.91860 -0.02126 -0.69548  0.117329  0.60753
## 24  0.03851  0.57991 -0.70662 -1.10691  0.533700 -0.10095
## 25  0.61401 -0.01017 -0.51875  0.25795 -1.076732 -0.24901
## 26  0.42275  0.07490  0.99168  0.37950 -0.249480  0.19511
## 27  0.54478 -0.12356  0.31406  0.41569  0.237094  0.09772
## 28  0.02044  0.43106 -0.71781  0.33274 -0.428303 -0.26278
## 29  1.23235  0.15293  0.10463 -0.07285  0.107338  0.07885
## 30  0.02214 -1.33353 -0.37271  0.38762 -0.062960 -0.48751
## 31 -0.71637  0.11204 -0.85950  0.24734 -0.066602 -0.17455
## 32  0.01756 -0.71899 -0.27770  0.01727  0.106473 -0.73172
## 33  0.78296 -0.75806  0.27650 -0.27767 -0.031297 -0.01895
## 34 -0.33590 -0.63554  1.07989 -0.07860 -0.526645  0.97933
## 35 -0.15761  0.49842  0.66426 -0.05091  0.331527  0.71522
## 36  0.53360  0.88375 -0.84613  0.03571 -0.630019 -0.08456
## 37  0.31300  0.02099 -0.05877 -0.57261  0.712218  0.05155
## 
## 
## Biplot scores for constraining variables
## 
##               dbRDA1   dbRDA2  dbRDA3   dbRDA4  dbRDA5   dbRDA6
## T_av         0.63982  0.25527  0.1662 -0.34674 -0.2289 -0.48193
## O2_sat_av    0.02371 -0.02938 -0.8854 -0.18557  0.2035 -0.13356
## Con_av      -0.25315 -0.27391 -0.1407 -0.76998 -0.4264  0.07657
## COD_av      -0.13235  0.15034  0.4807  0.08407  0.1430 -0.57341
## NH4._av     -0.40051  0.14159  0.7569 -0.21438 -0.2051 -0.37005
## Nt_av       -0.31714 -0.23687  0.2952 -0.33830 -0.5347 -0.40917
## pool_riffle -0.34609 -0.15029 -0.4384  0.27287 -0.5547 -0.42271
## meander     -0.33563  0.52703 -0.3800 -0.21165 -0.5422 -0.13159
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
##          Df SumOfSqs      F Pr(>F)
## Model     8   1.1664 1.2909   0.18
## Residual 28   3.1624
RsquareAdj(spe.rda)$adj.r.squared
## [1] 0.06072755
mod0 <- dbrda(spe.hel_bray ~ 1, env_select[,-c(9:10)])  # Model with intercept only  #edit_PH
mod1 <- dbrda(spe.hel_bray ~ ., env_select[,-c(9:10)])  # Model with all explanatory variables  #edit_PH
step.res <- ordiR2step(mod0, mod1, direction = "both",perm.max = 200)
## Step: R2.adj= 0 
## Call: spe.hel_bray ~ 1 
##  
##                  R2.adjusted
## <All variables>  0.060727547
## + T_av           0.036398792
## + NH4._av        0.020208612
## + meander        0.018502880
## + O2_sat_av      0.004277611
## + Con_av         0.001872668
## + pool_riffle    0.001742860
## <none>           0.000000000
## + Nt_av         -0.002060170
## + COD_av        -0.017968936
## 
##        Df    AIC      F Pr(>F)  
## + T_av  1 54.788 2.3599  0.052 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
step.res$anova  # Summary table
## NULL
plot(spe.rda, scaling = 1) # it is for technical reasons not possible to plot both site and species scores

summary(spe.rda)
## 
## Call:
## dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av +      NH4._av + Nt_av + pool_riffle + meander, data = env_select) 
## 
## Partitioning of squared Bray distance:
##               Inertia Proportion
## Total           4.329     1.0000
## Constrained     1.166     0.2695
## Unconstrained   3.162     0.7305
## 
## Eigenvalues, and their contribution to the squared Bray distance 
## 
## Importance of components:
##                       dbRDA1 dbRDA2  dbRDA3  dbRDA4   dbRDA5    dbRDA6
## Eigenvalue            0.5740 0.3372 0.17566 0.09415 0.042525 0.0020670
## Proportion Explained  0.1326 0.0779 0.04058 0.02175 0.009824 0.0004775
## Cumulative Proportion     NA     NA      NA      NA       NA        NA
##                         idbRDA1  idbRDA2   MDS1   MDS2   MDS3    MDS4    MDS5
## Eigenvalue            -0.014797 -0.04441 1.3249 0.8027 0.4680 0.31323 0.29121
## Proportion Explained   0.003418  0.01026 0.3061 0.1854 0.1081 0.07236 0.06727
## Cumulative Proportion        NA       NA     NA     NA     NA      NA      NA
##                          MDS6    MDS7    MDS8    MDS9   MDS10   MDS11   MDS12
## Eigenvalue            0.13644 0.11438 0.08987 0.07065 0.06425 0.02766 0.01576
## Proportion Explained  0.03152 0.02642 0.02076 0.01632 0.01484 0.00639 0.00364
## Cumulative Proportion      NA      NA      NA      NA      NA      NA      NA
##                          MDS13     MDS14      iMDS1     iMDS2    iMDS3
## Eigenvalue            0.011872 2.725e-04 -0.0020314 -0.007356 -0.01078
## Proportion Explained  0.002743 6.295e-05  0.0004693  0.001699  0.00249
## Cumulative Proportion       NA        NA         NA        NA       NA
##                           iMDS4    iMDS5     iMDS6     iMDS7    iMDS8     iMDS9
## Eigenvalue            -0.014744 -0.01935 -0.022495 -0.029849 -0.03342 -0.041489
## Proportion Explained   0.003406  0.00447  0.005197  0.006895  0.00772  0.009584
## Cumulative Proportion        NA       NA        NA        NA       NA        NA
##                         iMDS10   iMDS11   iMDS12   iMDS13   iMDS14
## Eigenvalue            -0.05984 -0.06439 -0.06995 -0.09116 -0.10203
## Proportion Explained   0.01382  0.01487  0.01616  0.02106  0.02357
## Cumulative Proportion       NA       NA       NA       NA       NA
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                       dbRDA1 dbRDA2 dbRDA3  dbRDA4  dbRDA5   dbRDA6  idbRDA1
## Eigenvalue            0.5740 0.3372 0.1757 0.09415 0.04252 0.002067 -0.01480
## Proportion Explained  0.4921 0.2891 0.1506 0.08072 0.03646 0.001772  0.01269
## Cumulative Proportion     NA     NA     NA      NA      NA       NA       NA
##                        idbRDA2
## Eigenvalue            -0.04441
## Proportion Explained   0.03807
## Cumulative Proportion       NA
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  3.533199 
## 
## 
## Site scores (weighted sums of species scores)
## 
##      dbRDA1   dbRDA2   dbRDA3    dbRDA4   dbRDA5   dbRDA6
## 1  -0.07530 -0.62580 -0.63602 -2.643196 -0.34427 -23.8946
## 2   0.87242  0.77832 -0.32068  2.074749 -0.25800  41.3754
## 3  -0.81094  1.29200 -0.57735 -0.708420 -2.35800  -5.0035
## 4  -0.97450  0.91051  1.79849 -2.192492 -1.99248  -4.7134
## 5  -0.25889  0.74559 -1.30995 -1.418227  0.33185 -21.9720
## 6   0.41640  1.01888 -0.34551 -2.318324  1.61924 -24.8949
## 7  -0.98286 -0.67495  0.62017  1.174218  1.55155  -0.8974
## 8   0.93807  0.62635 -3.10878 -0.765073  2.41981  -3.0764
## 9  -0.62816 -0.77245 -1.30540  1.806170  1.30543 -41.2679
## 10 -1.03727 -0.50089  0.62572  0.369727  1.00378   4.7819
## 11  1.10743 -0.99448 -0.66515 -0.071830  0.83402 -14.8528
## 12  0.61624  0.76828  0.87110  0.357836  0.12125  38.2552
## 13 -0.96088 -0.72453  0.47977  1.569869  1.43399  -6.8692
## 14  1.41128  2.03453 -1.67256  0.776255  0.27648  29.9435
## 15 -1.00448 -0.59302  0.65010  0.833950  1.44643   3.4902
## 16  1.26382  0.81623  3.19230  1.380128  1.12523 -10.4959
## 17 -0.92726 -0.73637  0.33774  1.892276  1.07787  -8.2400
## 18 -1.04523 -1.01280 -0.59052  2.897871  0.35870 -14.2865
## 19  0.89987 -1.54255 -0.54435  1.738188 -2.61939  32.8867
## 20 -0.14608 -0.34079  0.78630  1.316880  1.70249  -9.1246
## 21 -1.30124 -1.46318 -1.27370  3.739798  0.91928 -16.7232
## 22 -0.07485  1.73636  0.30728 -1.760677 -0.46041   6.1504
## 23 -0.18008 -0.64105  0.55428 -1.310793 -0.82437  -3.6579
## 24  0.53054  0.08088  1.31623 -3.529128 -1.96164  41.6937
## 25 -0.50901 -0.82222 -1.26480  2.167458 -6.71178 -22.4535
## 26  0.54874  0.92050  1.12409  0.009258  0.07714  34.0773
## 27  1.25804 -0.69768  2.37255 -1.697654 -3.24936 -16.2349
## 28 -0.86554  0.32492 -0.68347  0.704689 -0.16329  13.0277
## 29  1.82276  0.98120  0.33891 -1.852448  4.33581  -8.6587
## 30  0.72974 -0.59907 -1.00943 -0.201449  1.77124  -4.5186
## 31 -1.02904  0.19515  1.05102 -1.468951 -1.36610  23.1359
## 32 -0.13877 -0.06050 -0.19569 -1.899083 -1.10403 -31.4856
## 33  1.01209 -0.59332  2.58400 -0.318337 -1.43566  39.0635
## 34 -0.39319 -1.79743  0.20497 -0.362972 -0.59693 -19.8004
## 35 -0.89506  1.46497 -0.05622 -1.876722 -1.86120   3.5844
## 36  0.87386  1.17335 -3.50567  1.196922  1.51036  -9.4970
## 37 -0.06267 -0.67495 -0.14975  0.389535  2.08495  11.1529
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##      dbRDA1   dbRDA2   dbRDA3   dbRDA4    dbRDA5   dbRDA6
## 1  -0.16605 -0.69736 -0.45052  0.39949 -1.140870  0.26766
## 2  -0.07683 -0.33692  0.61638  0.50606 -0.381728  1.66676
## 3  -0.10071  0.76299  0.51759  0.02910  0.003366 -0.26259
## 4  -1.55548  0.90013  0.66933 -0.51241 -0.415318 -0.13093
## 5  -0.41914  0.27994 -0.66104 -1.10909  0.384267  0.50060
## 6   0.63090  0.30135 -0.07724  0.37780  0.857007 -0.33310
## 7  -0.15252 -0.25647 -0.20458 -0.38652  0.792931  0.48033
## 8  -0.14874  0.12938 -0.39866 -0.78891  0.083814  0.27677
## 9  -0.14359  0.11628  0.95439 -0.04629  0.609992 -0.64671
## 10 -0.74154 -0.65088  0.87072 -1.17151 -0.875610 -1.73263
## 11  0.77571 -0.48314  0.02321 -0.76003  0.259283 -0.01494
## 12 -0.44387 -0.28442  0.54938 -0.01813  1.062736 -0.71676
## 13 -0.25608 -0.61494 -0.47847  0.74269  0.180828 -0.14346
## 14 -0.23636  0.67472 -0.87185  0.72759  0.046993 -0.61306
## 15 -0.75508  0.19781  0.41035  0.89449 -0.744271 -0.25040
## 16  0.97946  1.50039  0.74911  0.12808 -0.657115 -0.30363
## 17 -1.19679 -0.18091 -0.90174  0.16093 -0.229189  0.70454
## 18  0.23197  0.39183 -0.42932 -0.82808 -0.242760  0.64991
## 19  0.80868 -0.55920 -0.22435  0.32426 -0.445610 -0.51746
## 20 -0.19206  0.12003  0.24728  0.99541  0.747110  0.40576
## 21 -0.38477  0.09371 -0.03632  1.07969  1.383306 -0.25819
## 22  0.02352  0.34014  0.07462  0.03662 -0.352806  0.35626
## 23  0.18716 -0.91860 -0.02126 -0.69548  0.117329  0.60753
## 24  0.03851  0.57991 -0.70662 -1.10691  0.533700 -0.10095
## 25  0.61401 -0.01017 -0.51875  0.25795 -1.076732 -0.24901
## 26  0.42275  0.07490  0.99168  0.37950 -0.249480  0.19511
## 27  0.54478 -0.12356  0.31406  0.41569  0.237094  0.09772
## 28  0.02044  0.43106 -0.71781  0.33274 -0.428303 -0.26278
## 29  1.23235  0.15293  0.10463 -0.07285  0.107338  0.07885
## 30  0.02214 -1.33353 -0.37271  0.38762 -0.062960 -0.48751
## 31 -0.71637  0.11204 -0.85950  0.24734 -0.066602 -0.17455
## 32  0.01756 -0.71899 -0.27770  0.01727  0.106473 -0.73172
## 33  0.78296 -0.75806  0.27650 -0.27767 -0.031297 -0.01895
## 34 -0.33590 -0.63554  1.07989 -0.07860 -0.526645  0.97933
## 35 -0.15761  0.49842  0.66426 -0.05091  0.331527  0.71522
## 36  0.53360  0.88375 -0.84613  0.03571 -0.630019 -0.08456
## 37  0.31300  0.02099 -0.05877 -0.57261  0.712218  0.05155
## 
## 
## Biplot scores for constraining variables
## 
##               dbRDA1   dbRDA2  dbRDA3   dbRDA4  dbRDA5   dbRDA6
## T_av         0.63982  0.25527  0.1662 -0.34674 -0.2289 -0.48193
## O2_sat_av    0.02371 -0.02938 -0.8854 -0.18557  0.2035 -0.13356
## Con_av      -0.25315 -0.27391 -0.1407 -0.76998 -0.4264  0.07657
## COD_av      -0.13235  0.15034  0.4807  0.08407  0.1430 -0.57341
## NH4._av     -0.40051  0.14159  0.7569 -0.21438 -0.2051 -0.37005
## Nt_av       -0.31714 -0.23687  0.2952 -0.33830 -0.5347 -0.40917
## pool_riffle -0.34609 -0.15029 -0.4384  0.27287 -0.5547 -0.42271
## meander     -0.33563  0.52703 -0.3800 -0.21165 -0.5422 -0.13159
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
##          Df SumOfSqs      F Pr(>F)
## Model     8   1.1664 1.2909  0.181
## Residual 28   3.1624
anova(spe.rda, by="term")
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
##             Df SumOfSqs      F Pr(>F)  
## T_av         1   0.2734 2.4210  0.061 .
## O2_sat_av    1   0.1377 1.2195  0.319  
## Con_av       1   0.1613 1.4283  0.246  
## COD_av       1   0.0676 0.5990  0.665  
## NH4._av      1   0.1709 1.5131  0.229  
## Nt_av        1   0.0501 0.4439  0.821  
## pool_riffle  1   0.0657 0.5818  0.692  
## meander      1   0.2396 2.1210  0.083 .
## Residual    28   3.1624                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova.cca(spe.rda, step=1000);
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
##          Df SumOfSqs      F Pr(>F)
## Model     8   1.1664 1.2909  0.153
## Residual 28   3.1624
anova.cca(spe.rda, step=1000, by="term");
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select)
##             Df SumOfSqs      F Pr(>F)  
## T_av         1   0.2734 2.4210  0.063 .
## O2_sat_av    1   0.1377 1.2195  0.287  
## Con_av       1   0.1613 1.4283  0.207  
## COD_av       1   0.0676 0.5990  0.677  
## NH4._av      1   0.1709 1.5131  0.212  
## Nt_av        1   0.0501 0.4439  0.792  
## pool_riffle  1   0.0657 0.5818  0.687  
## meander      1   0.2396 2.1210  0.087 .
## Residual    28   3.1624                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(spe.rda)$adj.r.squared;
## [1] 0.06072755
RsquareAdj(spe.rda)$r.squared
## [1] 0.2694548
# Generate a plot of the RDA
rownames(spe.rda$CCA$biplot) = c("temperature","O2","conductivity","COD","NH4","Nt","pool riffle pattern","meander")
ordiplot(spe.rda, choices=c(1,2), type="n", xlab="dbRDA 1 (49.21 %)", ylab="dbRDA 2 (28.91 %)")
points(spe.rda, "sites", pch=21, col="black", bg="white")
text(spe.rda, "cn", col="grey", cex=0.9, scaling=1)

# Save the plot to an object
figure4c = recordPlot()

8.1.2 Effect of space on component community structure

# Same for spatial predictors
spe.rda <- dbrda(spe.hel_bray ~ netcen + updist, data = env_select)
summary(spe.rda)
## 
## Call:
## dbrda(formula = spe.hel_bray ~ netcen + updist, data = env_select) 
## 
## Partitioning of squared Bray distance:
##               Inertia Proportion
## Total          4.3288     1.0000
## Constrained    0.5154     0.1191
## Unconstrained  3.8135     0.8809
## 
## Eigenvalues, and their contribution to the squared Bray distance 
## 
## Importance of components:
##                        dbRDA1  dbRDA2  MDS1   MDS2   MDS3  MDS4    MDS5    MDS6
## Eigenvalue            0.39474 0.12062 1.389 0.9772 0.7302 0.472 0.33207 0.19171
## Proportion Explained  0.09119 0.02787 0.321 0.2257 0.1687 0.109 0.07671 0.04429
## Cumulative Proportion      NA      NA    NA     NA     NA    NA      NA      NA
##                          MDS7   MDS8    MDS9   MDS10    MDS11    MDS12    MDS13
## Eigenvalue            0.14785 0.1195 0.09463 0.05735 0.032947 0.026192 0.015370
## Proportion Explained  0.03415 0.0276 0.02186 0.01325 0.007611 0.006051 0.003551
## Cumulative Proportion      NA     NA      NA      NA       NA       NA       NA
##                          MDS14     MDS15     MDS16      iMDS1      iMDS2
## Eigenvalue            0.011064 0.0023782 0.0016478 -0.0013922 -0.0028560
## Proportion Explained  0.002556 0.0005494 0.0003807  0.0003216  0.0006598
## Cumulative Proportion       NA        NA        NA         NA         NA
##                           iMDS3     iMDS4     iMDS5     iMDS6     iMDS7
## Eigenvalue            -0.005784 -0.014384 -0.016807 -0.017753 -0.021376
## Proportion Explained   0.001336  0.003323  0.003882  0.004101  0.004938
## Cumulative Proportion        NA        NA        NA        NA        NA
##                           iMDS8     iMDS9    iMDS10   iMDS11   iMDS12   iMDS13
## Eigenvalue            -0.027152 -0.034576 -0.036679 -0.04260 -0.04684 -0.05897
## Proportion Explained   0.006272  0.007988  0.008473  0.00984  0.01082  0.01362
## Cumulative Proportion        NA        NA        NA       NA       NA       NA
##                         iMDS14   iMDS15   iMDS16   iMDS17   iMDS18
## Eigenvalue            -0.06787 -0.07769 -0.08770 -0.10589 -0.12169
## Proportion Explained   0.01568  0.01795  0.02026  0.02446  0.02811
## Cumulative Proportion       NA       NA       NA       NA       NA
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                       dbRDA1 dbRDA2
## Eigenvalue            0.3947 0.1206
## Proportion Explained  0.7659 0.2341
## Cumulative Proportion 0.7659 1.0000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  3.533199 
## 
## 
## Site scores (weighted sums of species scores)
## 
##      dbRDA1   dbRDA2      MDS1      MDS2     MDS3     MDS4
## 1   0.03079  0.02086  0.443572 -0.509073 -0.07769 -1.05770
## 2   0.71832  0.43116  0.474971  0.635842 -0.29326  0.50678
## 3  -0.04721 -2.74094 -0.141670 -0.088012  1.14135 -0.12995
## 4   0.01772 -3.50459  0.115171 -1.313350  0.20198  0.84532
## 5  -0.32566 -1.42176  0.654893  0.016069  0.55155 -0.06800
## 6   0.72210 -1.64050  0.570204  0.088382  0.10581 -0.10658
## 7  -1.10558  0.09717 -0.369844 -0.272616 -0.16810  0.34656
## 8  -0.18673  0.52663  0.554834  1.248341  0.06126 -0.62322
## 9  -1.87557  1.15858 -0.639086  0.781404 -0.24620 -0.14224
## 10 -0.96878 -0.76984 -0.719617 -0.290543  0.49107  0.20709
## 11  0.11261  2.48593  0.221804  0.402044 -0.80474 -0.37386
## 12  1.04844 -0.21554  0.193892  0.092443  0.25970  0.86003
## 13 -1.24808  0.44331 -0.599846  0.028805 -0.19454  0.28743
## 14  1.25611 -0.45196  1.551121  0.523263  0.10082  0.21717
## 15 -1.01413 -0.32955 -0.980964  0.093579  0.09612 -0.04228
## 16  1.90373  0.67458  0.097707 -0.008746 -0.92598  1.56757
## 17 -1.37395  0.68301 -0.612465  0.141270 -0.23712  0.34682
## 18 -1.95385  1.26346 -0.581503  0.392340 -0.24748  0.26804
## 19  0.23442  3.46276 -0.164253  0.251482 -0.83790 -0.91794
## 20 -0.12218  0.38358 -0.532349  0.230205 -0.42232  0.20971
## 21 -2.52028  1.68572 -1.325649  1.074429  0.00633 -0.22347
## 22  1.09585 -2.60626  0.251395 -0.325519  0.60016  0.15487
## 23  0.13767 -0.04685 -0.057907 -0.573964 -0.42608 -0.44829
## 24  1.81459 -0.44404  0.571144 -1.006619 -0.20588 -0.86944
## 25 -1.70361  2.75887 -0.500433  0.366802  0.32907  0.38954
## 26  1.21825 -0.62815  0.148124  0.093074  0.05832  0.64989
## 27  2.04310  1.71402  0.001172 -0.871680 -0.94463 -0.70365
## 28 -0.62007 -1.14215 -0.376631  0.105274  1.18780 -0.07127
## 29  1.34657  0.13746  1.111295  0.401562 -0.91137  0.97202
## 30  0.26768  1.28419  0.333591  0.181043 -0.17537 -0.63533
## 31 -0.45809 -2.41232 -0.198776 -0.736562  0.25061  0.32623
## 32  0.11686 -0.75156  0.344145 -0.520102 -0.03041 -0.61279
## 33  1.57525  1.84844  0.082089 -0.494682 -0.76938  0.15659
## 34 -0.47482  1.79597 -0.639963 -0.480561 -0.17134 -0.79326
## 35  0.14388 -3.64825  0.208887 -0.768867  1.43955  0.43241
## 36  0.21397 -0.19845  0.700349  1.187955  1.06427 -0.60028
## 37 -0.01932  0.09703 -0.189405 -0.074710  0.14400 -0.32451
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##      dbRDA1   dbRDA2      MDS1      MDS2     MDS3     MDS4
## 1  -0.73416 -0.17667  0.443572 -0.509073 -0.07769 -1.05770
## 2   0.02228 -0.61344  0.474971  0.635842 -0.29326  0.50678
## 3   0.79019 -0.81424 -0.141670 -0.088012  1.14135 -0.12995
## 4  -0.64513 -0.18785  0.115171 -1.313350  0.20198  0.84532
## 5  -0.70290 -0.24187  0.654893  0.016069  0.55155 -0.06800
## 6   0.33049 -0.76784  0.570204  0.088382  0.10581 -0.10658
## 7  -1.13258  0.01401 -0.369844 -0.272616 -0.16810  0.34656
## 8   0.17308 -0.73123  0.554834  1.248341  0.06126 -0.62322
## 9   0.05980 -0.66540 -0.639086  0.781404 -0.24620 -0.14224
## 10  0.11807  0.66719 -0.719617 -0.290543  0.49107  0.20709
## 11  0.18170  0.62463  0.221804  0.402044 -0.80474 -0.37386
## 12  0.52300  0.62623  0.193892  0.092443  0.25970  0.86003
## 13 -0.67811 -0.25729 -0.599846  0.028805 -0.19454  0.28743
## 14 -0.93278 -0.10672  1.551121  0.523263  0.10082  0.21717
## 15  0.44549 -0.60608 -0.980964  0.093579  0.09612 -0.04228
## 16  0.90818 -0.07702  0.097707 -0.008746 -0.92598  1.56757
## 17 -0.70842 -0.20120 -0.612465  0.141270 -0.23712  0.34682
## 18 -1.15566  0.03171 -0.581503  0.392340 -0.24748  0.26804
## 19  0.32250  0.59444 -0.164253  0.251482 -0.83790 -0.91794
## 20  0.27139 -0.57771 -0.532349  0.230205 -0.42232  0.20971
## 21 -0.02181 -0.40054 -1.325649  1.074429  0.00633 -0.22347
## 22  0.78741 -0.94715  0.251395 -0.325519  0.60016  0.15487
## 23 -0.44652 -0.38594 -0.057907 -0.573964 -0.42608 -0.44829
## 24  0.29443 -0.76841  0.571144 -1.006619 -0.20588 -0.86944
## 25 -0.49871  1.02553 -0.500433  0.366802  0.32907  0.38954
## 26  0.62527 -0.36947  0.148124  0.093074  0.05832  0.64989
## 27  1.16557  0.23795  0.001172 -0.871680 -0.94463 -0.70365
## 28  0.46140  0.52227 -0.376631  0.105274  1.18780 -0.07127
## 29  0.06443  0.75038  1.111295  0.401562 -0.91137  0.97202
## 30 -0.05978  0.76030  0.333591  0.181043 -0.17537 -0.63533
## 31 -0.45764 -0.34376 -0.198776 -0.736562  0.25061  0.32623
## 32 -0.66474 -0.22103  0.344145 -0.520102 -0.03041 -0.61279
## 33  0.59703  0.43085  0.082089 -0.494682 -0.76938  0.15659
## 34  0.11174  0.84412 -0.639963 -0.480561 -0.17134 -0.79326
## 35  0.27483  0.73037  0.208887 -0.768867  1.43955  0.43241
## 36  0.26872  0.73193  0.700349  1.187955  1.06427 -0.60028
## 37  0.04198  0.86894 -0.189405 -0.074710  0.14400 -0.32451
## 
## 
## Biplot scores for constraining variables
## 
##         dbRDA1  dbRDA2 MDS1 MDS2 MDS3 MDS4
## netcen -0.8133  0.5819    0    0    0    0
## updist -0.8636 -0.5041    0    0    0    0
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ netcen + updist, data = env_select)
##          Df SumOfSqs      F Pr(>F)  
## Model     2   0.5154 2.2975  0.018 *
## Residual 34   3.8135                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(spe.rda)$adj.r.squared
## [1] 0.06723433
mod0 <- dbrda(spe.hel_bray ~ 1, env_select[,c(9:10)])  # Model with intercept only  #edit_PH
mod1 <- dbrda(spe.hel_bray ~ ., env_select[,c(9:10)])  # Model with all explanatory variables  #edit_PH
step.res <- ordiR2step(mod0, mod1, direction = "both",perm.max = 200)
## Step: R2.adj= 0 
## Call: spe.hel_bray ~ 1 
##  
##                 R2.adjusted
## <All variables>  0.06723433
## + updist         0.04867127
## + netcen         0.04317133
## <none>           0.00000000
## 
##          Df    AIC      F Pr(>F)  
## + updist  1 54.314 2.8418  0.022 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.04867127 
## Call: spe.hel_bray ~ updist 
##  
##                 R2.adjusted
## <All variables>  0.06723433
## + netcen         0.06723433
## <none>           0.04867127
## 
##          Df    AIC      F Pr(>F)
## + netcen  1 54.512 1.6965  0.172
step.res$anova  # Summary table
##                   R2.adj Df    AIC      F Pr(>F)  
## + updist        0.048671  1 54.314 2.8418  0.022 *
## <All variables> 0.067234                          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(spe.rda)$adj.r.squared
## [1] 0.06723433
anova.cca(spe.rda, step=1000);
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ netcen + updist, data = env_select)
##          Df SumOfSqs      F Pr(>F)  
## Model     2   0.5154 2.2975  0.028 *
## Residual 34   3.8135                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova.cca(spe.rda, step=1000, by="term");
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ netcen + updist, data = env_select)
##          Df SumOfSqs      F Pr(>F)  
## netcen    1   0.3019 2.6920  0.033 *
## updist    1   0.2134 1.9029  0.114  
## Residual 34   3.8135                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(spe.rda)$adj.r.squared;
## [1] 0.06723433
RsquareAdj(spe.rda)$r.squared
## [1] 0.1190546
# Generate a plot of the RDA
rownames(spe.rda$CCA$biplot) = c("network centrality","upstream distance")
ordiplot(spe.rda, choices=c(1,2), type="n", xlab="dbRDA 1 (76.59 %)", ylab="dbRDA 2 (23.41 %)")
points(spe.rda, "sites", pch=21, col="black", bg="white")
text(spe.rda, "cn", col="grey", cex=0.9, scaling=1, pos=4)
## Warning in arrows(0, 0, pts[, 1], pts[, 2], length = head.arrow, ...): "pos" is
## not a graphical parameter
## Warning in strwidth(labels, ...): "pos" is not a graphical parameter
## Warning in strheight(labels, ...): "pos" is not a graphical parameter

# Save the plot to an object
figure4d = recordPlot()

8.1.3 Variation partitioning

#Variation partitioning
spe.varpart1 <- varpart(spe.hel_bray, env_select[,1:8], env_select[,9:10])
plot(spe.varpart1,digits=2)

spe.varpart1
## 
## Partition of squared Bray distance in dbRDA 
## 
## Call: varpart(Y = spe.hel_bray, X = env_select[, 1:8], env_select[,
## 9:10])
## 
## Explanatory tables:
## X1:  env_select[, 1:8]
## X2:  env_select[, 9:10] 
## 
## No. of explanatory tables: 2 
## Total variation (SS): 4.3288 
## No. of observations: 37 
## 
## Partition table:
##                      Df R.squared Adj.R.squared Testable
## [a+c] = X1            8   0.26945       0.06073     TRUE
## [b+c] = X2            2   0.11905       0.06723     TRUE
## [a+b+c] = X1+X2      10   0.36716       0.12376     TRUE
## Individual fractions                                    
## [a] = X1|X2           8                 0.05653     TRUE
## [b] = X2|X1           2                 0.06303     TRUE
## [c]                   0                 0.00420    FALSE
## [d] = Residuals                         0.87624    FALSE
## ---
## Use function 'dbrda' to test significance of fractions of interest
anova.cca(dbrda(spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + Condition(netcen + updist),
                data=env_select), step=1000)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander + Condition(netcen + updist), data = env_select)
##          Df SumOfSqs      F Pr(>F)
## Model     8   1.0740 1.2742  0.181
## Residual 26   2.7395
anova.cca(dbrda(spe.hel_bray ~ netcen + updist+
                  Condition(T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander), data=env_select), step=1000)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = spe.hel_bray ~ netcen + updist + Condition(T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander), data = env_select)
##          Df SumOfSqs      F Pr(>F)  
## Model     2  0.42295 2.0071  0.069 .
## Residual 26  2.73945                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8.2 Infra-communities

# Infracommunities: Bray-Curtis dissimilarities are calculated at the individual host level Hellinger-transformed parasite data and then averaged within site
# A dummy parasite species is added to avoid problems with non-infected fishes
data_infra <- na.omit(data[,c(1,22:24,26:32)])
data_infra_disp <- dispweight(data_infra[,-1])
braycurtis <- vegdist(decostand(cbind(data_infra_disp,rep(1,nrow(data_infra))), na.rm=T, method="hellinger"), method="bray", na.rm=T)
meandist_bray <- meandist(braycurtis, data_infra[,1])

# Check whether Euclidean and Bray-Curtis distances are comparable
braycurtis <- vegdist(decostand(cbind(data_infra_disp,rep(1,nrow(data_infra))), na.rm=T, method="hellinger"), method="bray", na.rm=T)
meandist_bray <- meandist(braycurtis, data_infra[,1])
euc <- vegdist(decostand(cbind(data_infra_disp,rep(1,nrow(data_infra))), na.rm=T, method="hellinger"), method="euc", na.rm=T)
meandist_euc <- meandist(euc, data_infra[,1])
plot(meandist_bray[1:37,1:37], meandist_euc[1:37,1:37])

mantel(meandist_bray[1:37,1:37], meandist_euc[1:37,1:37])
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = meandist_bray[1:37, 1:37], ydis = meandist_euc[1:37,      1:37]) 
## 
## Mantel statistic r: 0.9906 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.177 0.234 0.297 0.346 
## Permutation: free
## Number of permutations: 999
model_adonis = adonis(meandist_bray ~ avlength + avcondition + T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + netcen + 
                       updist, data=environment2)
## 'adonis' will be deprecated: use 'adonis2' instead
model_adonis$aov.tab
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##             Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)   
## avlength     1   0.05123 0.051226  6.1160 0.11570  0.009 **
## avcondition  1   0.00393 0.003927  0.4688 0.00887  0.653   
## T_av         1   0.00874 0.008741  1.0437 0.01974  0.393   
## O2_sat_av    1   0.02116 0.021158  2.5262 0.04779  0.121   
## Con_av       1   0.05861 0.058606  6.9971 0.13237  0.004 **
## COD_av       1   0.03543 0.035428  4.2298 0.08002  0.026 * 
## NH4._av      1   0.00275 0.002754  0.3288 0.00622  0.708   
## Nt_av        1   0.01049 0.010488  1.2521 0.02369  0.328   
## pool_riffle  1   0.00491 0.004914  0.5867 0.01110  0.595   
## meander      1   0.00959 0.009588  1.1447 0.02166  0.335   
## netcen       1   0.02319 0.023187  2.7684 0.05237  0.084 . 
## updist       1   0.01171 0.011707  1.3977 0.02644  0.278   
## Residuals   24   0.20102 0.008376         0.45403          
## Total       36   0.44274                  1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# environmental variables
env_select <- environment2[,c("T_av", "O2_sat_av", "Con_av", "COD_av", "NH4._av", "Nt_av", "pool_riffle", "meander", "netcen", "updist")]
env_select$pool_riffle <- as.numeric(env_select$pool_riffle)
env_select$meander <- as.numeric(env_select$meander)

pca <- prcomp(env_select, scale.=T)
summary(pca)
## Importance of components:
##                           PC1    PC2    PC3    PC4     PC5     PC6     PC7
## Standard deviation     1.7124 1.5545 1.1221 1.0140 0.88807 0.79463 0.56647
## Proportion of Variance 0.2933 0.2416 0.1259 0.1028 0.07887 0.06314 0.03209
## Cumulative Proportion  0.2933 0.5349 0.6608 0.7636 0.84248 0.90563 0.93771
##                            PC8     PC9    PC10
## Standard deviation     0.50483 0.46939 0.38429
## Proportion of Variance 0.02549 0.02203 0.01477
## Cumulative Proportion  0.96320 0.98523 1.00000
plot(pca)

biplot(pca)

8.2.1 Effect of environment on infracommunity structure

# Assess the effect of environmental variables on parasite infracommunity dissimilarities using distance based RDA
spe.rda <- dbrda(meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander, data = env_select, tidy=TRUE)
summary(spe.rda)
## 
## Call:
## dbrda(formula = meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av +      NH4._av + Nt_av + pool_riffle + meander, data = env_select,      tidy = TRUE) 
## 
## Partitioning of squared Unknown distance:
##               Inertia Proportion
## Total          1.5898     1.0000
## Constrained    0.4504     0.2833
## Unconstrained  1.1395     0.7167
## 
## Eigenvalues, and their contribution to the squared Unknown distance 
## 
## Importance of components:
##                       dbRDA1  dbRDA2  dbRDA3  dbRDA4  dbRDA5  dbRDA6  dbRDA7
## Eigenvalue            0.2196 0.05240 0.04219 0.03708 0.02865 0.02593 0.02365
## Proportion Explained  0.1381 0.03296 0.02654 0.02333 0.01802 0.01631 0.01487
## Cumulative Proportion 0.1381 0.17107 0.19761 0.22093 0.23895 0.25527 0.27014
##                        dbRDA8   MDS1    MDS2    MDS3    MDS4    MDS5    MDS6
## Eigenvalue            0.02089 0.1762 0.13809 0.10605 0.06557 0.06183 0.05722
## Proportion Explained  0.01314 0.1108 0.08686 0.06671 0.04125 0.03889 0.03599
## Cumulative Proportion 0.28328 0.3941 0.48098 0.54769 0.58894 0.62783 0.66382
##                          MDS7    MDS8    MDS9   MDS10   MDS11   MDS12   MDS13
## Eigenvalue            0.04417 0.04047 0.03937 0.03784 0.03506 0.03412 0.03273
## Proportion Explained  0.02778 0.02545 0.02476 0.02380 0.02205 0.02146 0.02058
## Cumulative Proportion 0.69160 0.71706 0.74182 0.76562 0.78767 0.80913 0.82972
##                         MDS14   MDS15   MDS16   MDS17   MDS18   MDS19   MDS20
## Eigenvalue            0.02795 0.02650 0.02630 0.02513 0.02348 0.02044 0.01991
## Proportion Explained  0.01758 0.01667 0.01654 0.01580 0.01477 0.01286 0.01252
## Cumulative Proportion 0.84730 0.86397 0.88051 0.89632 0.91109 0.92394 0.93647
##                         MDS21   MDS22   MDS23    MDS24    MDS25    MDS26
## Eigenvalue            0.01975 0.01836 0.01714 0.015343 0.013426 0.008684
## Proportion Explained  0.01242 0.01155 0.01078 0.009651 0.008445 0.005462
## Cumulative Proportion 0.94889 0.96044 0.97122 0.980869 0.989313 0.994775
##                          MDS27     MDS28
## Eigenvalue            0.007381 0.0009256
## Proportion Explained  0.004642 0.0005822
## Cumulative Proportion 0.999418 1.0000000
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                       dbRDA1 dbRDA2  dbRDA3  dbRDA4  dbRDA5  dbRDA6  dbRDA7
## Eigenvalue            0.2196 0.0524 0.04219 0.03708 0.02865 0.02593 0.02365
## Proportion Explained  0.4875 0.1164 0.09368 0.08234 0.06361 0.05759 0.05251
## Cumulative Proportion 0.4875 0.6039 0.69757 0.77991 0.84353 0.90111 0.95362
##                        dbRDA8
## Eigenvalue            0.02089
## Proportion Explained  0.04638
## Cumulative Proportion 1.00000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  2.750505 
## 
## 
## Site scores (weighted sums of species scores)
## 
##            dbRDA1    dbRDA2   dbRDA3    dbRDA4   dbRDA5    dbRDA6
## SITE 1  -0.193401  0.063457 -0.61764 -0.312982  0.78782 -0.875368
## SITE 11 -0.424501 -0.396544 -0.49524 -0.201551  1.04886  0.509773
## SITE 12 -0.702700  0.085074 -0.27366  0.063060  0.10379 -0.004917
## SITE 13 -0.007236 -0.112556 -0.10006  0.582155  0.14248  0.722762
## SITE 14 -0.125146  0.085307 -0.29498  0.469582 -0.32057  0.196422
## SITE 15 -0.265297 -0.146273  0.10540  0.221087 -0.44493  0.501526
## SITE 16  0.017380  0.097078 -0.34338  0.027222 -0.23381  0.253728
## SITE 17 -0.143898  0.261154 -0.74786  0.425484 -1.19344 -1.632875
## SITE 18  0.694672  0.250138 -0.19428  0.116751 -0.02255  0.197381
## SITE 19  2.013583  0.700678 -0.26310  0.241350  0.28881 -0.781458
## SITE 2  -0.269160 -0.474927 -0.05506 -0.658268 -0.75242 -0.133679
## SITE 20  0.184271 -1.120553  0.05670  0.218183 -0.85772  0.307942
## SITE 21 -0.478165  0.259320 -0.31414 -0.196556  0.20691  0.132081
## SITE 22 -0.680264  0.661630  0.17236  0.300656  0.03888  0.186979
## SITE 23  0.107511  0.117471 -0.68571  0.024261  0.82469 -0.140649
## SITE 24 -0.509800 -0.585796  1.11284  0.791260  0.27686 -0.880620
## SITE 26  0.023452  0.520302 -0.75929 -0.348237  0.26634 -0.019943
## SITE 28  0.081124  0.364618 -0.32704 -0.131404 -0.21902 -0.874784
## SITE 29  0.069397  0.046767  0.74532 -1.644458  0.22297  0.183059
## SITE 3  -0.655376 -0.078746 -0.30374  0.077802  0.17383  0.171647
## SITE 30 -0.750306  0.186160 -0.53020 -0.087275  0.13772  0.020057
## SITE 31 -0.003379 -0.217568  0.16729  0.877596  0.12615 -0.182314
## SITE 32  0.525097 -0.234069 -0.10169  0.176823 -0.13431  0.273050
## SITE 33 -0.376746 -0.057152 -0.09869  0.386971 -0.79147 -0.032496
## SITE 34 -0.427146  0.315945 -0.34228 -0.588588  0.30660 -0.646267
## SITE 35  0.144043 -0.736009  1.21209  0.693519  0.92859  0.065714
## SITE 36 -0.221245 -1.014780  0.73058 -0.257720  0.04394  0.325929
## SITE 38  0.260121  0.548181 -0.12839  0.009881  0.12554  0.030680
## SITE 39 -0.912898  0.029942 -0.24674 -0.094014  0.07078  0.062174
## SITE 4  -0.039914 -0.160458 -0.52023 -1.243446  0.02703  0.155617
## SITE 40  0.437382  0.347286 -0.44241  0.551293  0.15660  0.371577
## SITE 41  1.061747  0.460820  0.41788 -0.329177 -0.20931  1.117435
## SITE 42  0.492754 -2.142555  1.38058 -0.595630 -0.70424 -0.924924
## SITE 5   0.486220 -0.097372 -0.04666  0.145692  0.78718  0.142772
## SITE 6   0.366097  0.007433  0.29450  1.020707  0.06306  1.202313
## SITE 7  -0.109776  2.157876  1.76424 -0.745567 -0.40355 -0.618902
## SITE 9   0.331503  0.008719  0.07268  0.013537 -0.86805  0.618579
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##            dbRDA1   dbRDA2    dbRDA3   dbRDA4   dbRDA5    dbRDA6
## SITE 1  -0.013706  0.06161 -0.566299 -0.56263  0.68337 -0.799399
## SITE 11 -0.498049 -0.62259 -0.387895 -0.14443  0.86286  0.337683
## SITE 12 -0.449977 -0.41292 -0.381703  0.39230 -0.03446 -0.494715
## SITE 13  0.646634  0.78752 -0.170195  0.73514  0.33780  0.894640
## SITE 14 -0.022889  0.21915 -0.559986  0.43068 -0.75852 -0.249161
## SITE 15 -0.477115 -0.18072  0.381861 -0.03891 -0.43609  0.261592
## SITE 16 -0.003779 -0.10419 -0.224779  0.04360 -0.47874  0.521790
## SITE 17  0.014080 -0.01398 -0.393448  0.37555 -0.45946 -0.612432
## SITE 18  0.433470 -0.38610  0.311804  0.49069 -0.10431  0.385204
## SITE 19  1.797552 -0.01477 -0.068343  0.14213 -0.02578 -0.648483
## SITE 2   0.097800 -0.47743  0.114446 -0.54828 -0.66048 -0.040147
## SITE 20  0.233513 -0.56819 -0.452230  0.18420 -0.43596  0.342081
## SITE 21  0.115261  0.39505  0.012004 -0.45991  0.25232  0.602002
## SITE 22 -0.271628  0.94479  0.262170  0.24479  0.06794 -0.038339
## SITE 23 -0.268975 -0.01406 -0.687107 -0.14113  0.84944 -0.150595
## SITE 24 -0.470197 -0.06626  1.074097  0.59004  0.26654 -0.627883
## SITE 26 -0.127347  0.66408 -0.992149 -0.04145  0.27143  0.182401
## SITE 28 -0.129789  0.14568  0.024292  0.34127 -0.32756 -0.565381
## SITE 29  0.193897  0.19946  0.677019 -0.93440  0.11201  0.170499
## SITE 3  -0.678081 -0.39312 -0.314876  0.16522  0.26479  0.467313
## SITE 30 -0.606829 -0.31255 -0.451988  0.29629 -0.13401  0.461075
## SITE 31 -0.139331 -0.12436  0.009754  0.63248  0.30366 -0.804481
## SITE 32  0.410116 -0.40478 -0.132639 -0.21685 -0.23921 -0.040507
## SITE 33 -0.074343  0.29536 -0.134441  0.39097 -1.04116 -0.336434
## SITE 34 -0.226204  0.29969  0.198989 -0.65741  0.28977 -0.726017
## SITE 35  0.238182 -0.36205  0.851383  0.44887  0.64367  0.065366
## SITE 36 -0.212694 -0.42956  0.325410 -0.02991  0.12247 -0.002134
## SITE 38 -0.290761  0.68825  0.099402 -0.20120  0.08529 -0.165182
## SITE 39 -0.424749 -0.34091  0.653590 -0.45291 -0.31773  0.026280
## SITE 4   0.249290 -0.18203 -0.478520 -1.12533  0.02918  0.073437
## SITE 40  0.011975  0.77178 -0.349485  0.10183  0.04876 -0.016774
## SITE 41  0.562714  0.30100  0.218979 -0.62325 -0.17919  0.463349
## SITE 42  0.137879 -0.72177  0.095094 -0.49959 -0.18024 -0.377517
## SITE 5   0.741324 -0.45934  0.254706  0.55877  0.92073  0.108541
## SITE 6  -0.173401 -0.11045  0.214506  0.34726  0.12331  0.892139
## SITE 7  -0.341411  1.02435  0.770083 -0.17660 -0.04056 -0.023937
## SITE 9   0.017570 -0.09563  0.196492 -0.05789 -0.68189  0.464125
## 
## 
## Biplot scores for constraining variables
## 
##               dbRDA1   dbRDA2   dbRDA3   dbRDA4   dbRDA5   dbRDA6
## T_av         0.19654  0.04389  0.79605 -0.13924 -0.28696 -0.21269
## O2_sat_av   -0.24302  0.51209 -0.29524 -0.44557 -0.62386 -0.02169
## Con_av       0.69249  0.28619 -0.13247 -0.04714 -0.14852 -0.19236
## COD_av       0.39564 -0.11996  0.29825  0.49441  0.08302 -0.03404
## NH4._av      0.62169 -0.12508  0.08817  0.43757  0.26414  0.12573
## Nt_av        0.61125 -0.07446 -0.30978 -0.09649  0.17024 -0.41824
## pool_riffle  0.25108  0.73042 -0.13229 -0.45593  0.34283 -0.02248
## meander     -0.04274  0.63037 -0.22360  0.31145  0.06195 -0.62153
rownames(spe.rda$CCA$biplot) = c("temperature","O2","conductivity","COD","NH4","Nt","pool riffle pattern","meander")

# Generate a plot of the RDA
ordiplot(spe.rda, choices=c(1,2), type="n", xlab="dbRDA 1 (48.75 %)", ylab="dbRDA 2 (11.64 %)")
points(spe.rda, "sites", pch=21, col="black", bg="white")
text(spe.rda, "cn", col="grey", cex=0.9, scaling=1)

# Save the plot to an object
figure4a = recordPlot()

anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select, tidy = TRUE)
##          Df SumOfSqs      F Pr(>F)   
## Model     8  0.45036 1.3833  0.007 **
## Residual 28  1.13946                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(spe.rda, by="term")
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander, data = env_select, tidy = TRUE)
##             Df SumOfSqs      F Pr(>F)    
## T_av         1  0.04338 1.0661  0.308    
## O2_sat_av    1  0.04859 1.1940  0.237    
## Con_av       1  0.12820 3.1503  0.001 ***
## COD_av       1  0.06629 1.6290  0.056 .  
## NH4._av      1  0.02849 0.7002  0.835    
## Nt_av        1  0.03286 0.8076  0.688    
## pool_riffle  1  0.04004 0.9838  0.397    
## meander      1  0.06250 1.5358  0.057 .  
## Residual    28  1.13946                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RsquareAdj(spe.rda)$adj.r.squared
## [1] 0.07850103
mod0 <- dbrda(meandist_bray ~ 1, env_select)  # Model with intercept only  #edit_PH
mod1 <- dbrda(meandist_bray ~ ., env_select)  # Model with all explanatory variables  #edit_PH
step.res <- ordiR2step(mod0, mod1, direction = "both",perm.max = 200)
## Step: R2.adj= 0 
## Call: meandist_bray ~ 1 
##  
##                   R2.adjusted
## <All variables>  0.0898822750
## + Con_av         0.0492715387
## + NH4._av        0.0376899022
## + Nt_av          0.0353578268
## + COD_av         0.0097867139
## + updist         0.0092890919
## + pool_riffle    0.0070398050
## + netcen         0.0034499960
## + O2_sat_av      0.0031240320
## <none>           0.0000000000
## + T_av          -0.0005031246
## + meander       -0.0037253892
## 
##          Df    AIC      F Pr(>F)   
## + Con_av  1 17.228 2.8657  0.004 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.04927154 
## Call: meandist_bray ~ Con_av 
##  
##                 R2.adjusted
## <All variables>  0.08988227
## + COD_av         0.08426096
## + NH4._av        0.07058808
## + updist         0.06800179
## + O2_sat_av      0.06266812
## + netcen         0.05455352
## + meander        0.05428744
## + Nt_av          0.05273614
## <none>           0.04927154
## + pool_riffle    0.04923534
## + T_av           0.04694754
## 
##          Df    AIC      F Pr(>F)   
## + COD_av  1 16.768 2.3373  0.004 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Step: R2.adj= 0.08426096 
## Call: meandist_bray ~ Con_av + COD_av 
##  
##                 R2.adjusted
## + updist         0.09961037
## + meander        0.09248673
## <All variables>  0.08988227
## + netcen         0.08626581
## + pool_riffle    0.08559576
## <none>           0.08426096
## + T_av           0.08115910
## + O2_sat_av      0.08092067
## + Nt_av          0.07979657
## + NH4._av        0.07591794
step.res$anova  # Summary table
##                   R2.adj Df    AIC      F Pr(>F)   
## + Con_av        0.049272  1 17.228 2.8657  0.004 **
## + COD_av        0.084261  1 16.768 2.3373  0.004 **
## <All variables> 0.089882                           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
spe.rda <- dbrda(meandist_bray ~  Con_av + COD_av, env_select)
plot(spe.rda, scaling = 1) # it is for technical reasons not possible to plot both site and species scores

summary(spe.rda)
## 
## Call:
## dbrda(formula = meandist_bray ~ Con_av + COD_av, data = env_select) 
## 
## Partitioning of squared Unknown distance:
##               Inertia Proportion
## Total          1.5898     1.0000
## Constrained    0.2148     0.1351
## Unconstrained  1.3750     0.8649
## 
## Eigenvalues, and their contribution to the squared Unknown distance 
## 
## Importance of components:
##                       dbRDA1  dbRDA2   MDS1    MDS2   MDS3    MDS4    MDS5
## Eigenvalue            0.1863 0.02858 0.1972 0.15873 0.1207 0.08661 0.07076
## Proportion Explained  0.1172 0.01798 0.1240 0.09984 0.0759 0.05447 0.04451
## Cumulative Proportion     NA      NA     NA      NA     NA      NA      NA
##                          MDS6    MDS7    MDS8    MDS9   MDS10   MDS11   MDS12
## Eigenvalue            0.06231 0.04802 0.04613 0.04387 0.04231 0.03826 0.03703
## Proportion Explained  0.03919 0.03021 0.02901 0.02760 0.02661 0.02406 0.02329
## Cumulative Proportion      NA      NA      NA      NA      NA      NA      NA
##                         MDS13   MDS14   MDS15   MDS16   MDS17   MDS18   MDS19
## Eigenvalue            0.03557 0.03355 0.03200 0.02808 0.02797 0.02575 0.02487
## Proportion Explained  0.02237 0.02110 0.02013 0.01766 0.01759 0.01620 0.01564
## Cumulative Proportion      NA      NA      NA      NA      NA      NA      NA
##                         MDS20   MDS21   MDS22   MDS23   MDS24   MDS25   MDS26
## Eigenvalue            0.02377 0.02238 0.02182 0.02104 0.01990 0.01880 0.01782
## Proportion Explained  0.01495 0.01408 0.01372 0.01323 0.01252 0.01183 0.01121
## Cumulative Proportion      NA      NA      NA      NA      NA      NA      NA
##                         MDS27   MDS28    MDS29    MDS30    MDS31    MDS32
## Eigenvalue            0.01719 0.01622 0.011778 0.010499 0.007903 0.006958
## Proportion Explained  0.01081 0.01020 0.007408 0.006604 0.004971 0.004376
## Cumulative Proportion      NA      NA       NA       NA       NA       NA
##                          MDS33     iMDS1
## Eigenvalue            0.005670 -0.006426
## Proportion Explained  0.003566  0.004042
## Cumulative Proportion       NA        NA
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                       dbRDA1  dbRDA2
## Eigenvalue            0.1863 0.02858
## Proportion Explained  0.8670 0.13304
## Cumulative Proportion 0.8670 1.00000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  2.750505 
## 
## 
## Site scores (weighted sums of species scores)
## 
##            dbRDA1   dbRDA2     MDS1     MDS2     MDS3      MDS4
## SITE 1  -0.215744  0.44381  0.01302  0.31819 -0.28116 -0.178171
## SITE 11 -0.507248  0.82144 -0.09719 -0.28198 -0.04679  0.342258
## SITE 12 -0.746029  0.15775  0.31672  0.25334 -0.08234 -0.151368
## SITE 13 -0.005847 -0.21986 -0.02027  0.91847 -0.12233 -0.332031
## SITE 14 -0.101987  0.28502  0.01314  0.55873 -0.09066 -0.223950
## SITE 15 -0.290401 -0.51714 -0.23614 -0.10408 -0.36670 -0.263022
## SITE 16  0.016467  0.22318  0.14693 -0.09879 -0.21010  0.156923
## SITE 17 -0.053594  0.78979  0.68552  0.22623  0.55763  0.068555
## SITE 18  0.741297 -0.34140  0.34156 -0.78781 -0.33039  0.372025
## SITE 19  2.204147  0.31903  0.07576 -0.77176 -0.32619  0.237586
## SITE 2  -0.334822  0.53901  0.20181  0.23760  0.25255  0.460075
## SITE 20  0.132501 -1.43824 -0.36307 -0.15683  0.91526  0.652586
## SITE 21 -0.518164  0.15765  0.51672  0.22690 -0.05180  0.071838
## SITE 22 -0.681351 -0.71637  0.44536  0.61813  0.24285 -0.440713
## SITE 23  0.085803 -0.07384 -0.18742 -0.58764 -0.50280  0.146049
## SITE 24 -0.506368 -1.41199 -0.02580  0.53865  0.25884  0.040947
## SITE 26  0.017028  0.61108  0.37526 -0.36715 -0.10492  0.425754
## SITE 28  0.121498  0.47609  0.59944 -0.22080 -0.07107  0.465237
## SITE 29 -0.014961  0.35365  0.04914 -0.10076 -0.08378  0.641604
## SITE 3  -0.711132 -0.09645  0.08469 -0.08718 -0.21517  0.038217
## SITE 30 -0.812113 -0.12175  0.55582 -0.16918 -0.25866  0.087617
## SITE 31  0.043616 -0.61571 -0.31459  0.47663 -0.16031 -0.651535
## SITE 32  0.562390  0.35603 -0.38737  0.06837 -0.20621 -0.072343
## SITE 33 -0.372173  0.03501  0.10319  0.85337 -0.01445 -0.328696
## SITE 34 -0.466651  0.48378  0.50738  0.02349 -0.08407  0.306511
## SITE 35  0.163917 -1.66724 -0.04520  0.22652  0.57157  0.198776
## SITE 36 -0.296043 -0.76361 -0.59812  0.40106 -0.04152 -0.249681
## SITE 38  0.302702  0.34509 -0.16274 -0.49724 -0.25040 -0.005416
## SITE 39 -0.978171  0.31142  0.36167  0.40699  0.02355  0.020399
## SITE 4  -0.141634  0.44343 -0.02774 -0.11080 -0.07834  0.180628
## SITE 40  0.496377 -0.01681 -0.38350 -0.04028 -0.50240 -0.167352
## SITE 41  1.117637  0.21808 -0.82799 -0.60061 -0.49878 -0.637014
## SITE 42  0.437517 -0.97421 -1.46550  0.08641  1.33505  0.880493
## SITE 5   0.522987  0.17736  0.27513  0.18757 -0.04206  0.152978
## SITE 6   0.423826 -0.07352 -0.91744 -0.10155 -0.42928 -0.789691
## SITE 7   0.013170  1.17661  0.57227 -1.15454  1.52741 -1.616432
## SITE 9   0.351552  0.32386 -0.18044 -0.38769 -0.23204  0.160360
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##            dbRDA1   dbRDA2     MDS1     MDS2     MDS3      MDS4
## SITE 1  -0.079483  0.44275  0.01302  0.31819 -0.28116 -0.178171
## SITE 11 -0.692558  0.69501 -0.09719 -0.28198 -0.04679  0.342258
## SITE 12 -0.387516 -0.29270  0.31672  0.25334 -0.08234 -0.151368
## SITE 13  0.692033  0.28074 -0.02027  0.91847 -0.12233 -0.332031
## SITE 14  0.273072  0.58517  0.01314  0.55873 -0.09066 -0.223950
## SITE 15 -0.503973 -0.49304 -0.23614 -0.10408 -0.36670 -0.263022
## SITE 16 -0.041630  0.31054  0.14693 -0.09879 -0.21010  0.156923
## SITE 17  0.286975  0.22462  0.68552  0.22623  0.55763  0.068555
## SITE 18  0.355899 -0.83587  0.34156 -0.78781 -0.33039  0.372025
## SITE 19  1.702755 -0.40556  0.07576 -0.77176 -0.32619  0.237586
## SITE 2  -0.068404  0.54101  0.20181  0.23760  0.25255  0.460075
## SITE 20  0.005618 -0.62143 -0.36307 -0.15683  0.91526  0.652586
## SITE 21 -0.085484 -0.08616  0.51672  0.22690 -0.05180  0.071838
## SITE 22 -0.015359 -0.63175  0.44536  0.61813  0.24285 -0.440713
## SITE 23 -0.475351 -0.06131 -0.18742 -0.58764 -0.50280  0.146049
## SITE 24 -0.138656 -0.63048 -0.02580  0.53865  0.25884  0.040947
## SITE 26 -0.086478  0.70999  0.37526 -0.36715 -0.10492  0.425754
## SITE 28  0.227529  0.47771  0.59944 -0.22080 -0.07107  0.465237
## SITE 29 -0.053389  0.16940  0.04914 -0.10076 -0.08378  0.641604
## SITE 3  -0.805011 -0.44279  0.08469 -0.08718 -0.21517  0.038217
## SITE 30 -0.711357 -0.97402  0.55582 -0.16918 -0.25866  0.087617
## SITE 31  0.197349 -0.41392 -0.31459  0.47663 -0.16031 -0.651535
## SITE 32  0.342870  0.50733 -0.38737  0.06837 -0.20621 -0.072343
## SITE 33  0.260419  0.25833  0.10319  0.85337 -0.01445 -0.328696
## SITE 34 -0.213627  0.34583  0.50738  0.02349 -0.08407  0.306511
## SITE 35  0.299656 -0.71657 -0.04520  0.22652  0.57157  0.198776
## SITE 36 -0.265626 -0.39128 -0.59812  0.40106 -0.04152 -0.249681
## SITE 38 -0.177490  0.11612 -0.16274 -0.49724 -0.25040 -0.005416
## SITE 39 -0.502932  0.21901  0.36167  0.40699  0.02355  0.020399
## SITE 4  -0.213933  0.31087 -0.02774 -0.11080 -0.07834  0.180628
## SITE 40  0.157625 -0.02847 -0.38350 -0.04028 -0.50240 -0.167352
## SITE 41  0.308647  0.03805 -0.82799 -0.60061 -0.49878 -0.637014
## SITE 42 -0.043511  0.13560 -1.46550  0.08641  1.33505  0.880493
## SITE 5   0.774214 -0.17770  0.27513  0.18757 -0.04206  0.152978
## SITE 6  -0.221478  0.23059 -0.91744 -0.10155 -0.42928 -0.789691
## SITE 7  -0.080948  0.35692  0.57227 -1.15454  1.52741 -1.616432
## SITE 9  -0.020466  0.24745 -0.18044 -0.38769 -0.23204  0.160360
## 
## 
## Biplot scores for constraining variables
## 
##        dbRDA1  dbRDA2 MDS1 MDS2 MDS3 MDS4
## Con_av 0.7628  0.6467    0    0    0    0
## COD_av 0.4413 -0.8974    0    0    0    0
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ Con_av + COD_av, data = env_select)
##          Df SumOfSqs      F Pr(>F)    
## Model     2  0.21484 2.6563  0.001 ***
## Residual 34  1.37498                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(spe.rda, by="term")
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ Con_av + COD_av, data = env_select)
##          Df SumOfSqs      F Pr(>F)    
## Con_av    1  0.12032 2.9752  0.001 ***
## COD_av    1  0.09452 2.3373  0.011 *  
## Residual 34  1.37498                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8.2.2 Effect of space on infracommunity structure

spe.rda <- dbrda(meandist_bray ~ netcen + updist, data = env_select)
anova(spe.rda)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ netcen + updist, data = env_select)
##          Df SumOfSqs      F Pr(>F)
## Model     2  0.10742 1.2319  0.127
## Residual 34  1.48240
RsquareAdj(spe.rda)$adj.r.squared
## [1] 0.01271734
anova.cca(spe.rda, step=1000);
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ netcen + updist, data = env_select)
##          Df SumOfSqs      F Pr(>F)
## Model     2  0.10742 1.2319  0.121
## Residual 34  1.48240
anova.cca(spe.rda, step=1000, by="term");
## Permutation test for dbrda under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ netcen + updist, data = env_select)
##          Df SumOfSqs      F Pr(>F)
## netcen    1  0.04949 1.1352  0.230
## updist    1  0.05792 1.3285  0.132
## Residual 34  1.48240
RsquareAdj(spe.rda)$adj.r.squared;
## [1] 0.01271734
RsquareAdj(spe.rda)$r.squared
## [1] 0.06756637
summary(spe.rda)
## 
## Call:
## dbrda(formula = meandist_bray ~ netcen + updist, data = env_select) 
## 
## Partitioning of squared Unknown distance:
##               Inertia Proportion
## Total          1.5898    1.00000
## Constrained    0.1074    0.06757
## Unconstrained  1.4824    0.93243
## 
## Eigenvalues, and their contribution to the squared Unknown distance 
## 
## Importance of components:
##                        dbRDA1  dbRDA2   MDS1   MDS2    MDS3    MDS4    MDS5
## Eigenvalue            0.05871 0.04871 0.3207 0.1666 0.10533 0.08746 0.06478
## Proportion Explained  0.03693 0.03064 0.2017 0.1048 0.06625 0.05501 0.04075
## Cumulative Proportion      NA      NA     NA     NA      NA      NA      NA
##                          MDS6    MDS7    MDS8    MDS9   MDS10   MDS11   MDS12
## Eigenvalue            0.06249 0.04788 0.04637 0.04324 0.04235 0.03756 0.03592
## Proportion Explained  0.03931 0.03012 0.02917 0.02720 0.02664 0.02362 0.02260
## Cumulative Proportion      NA      NA      NA      NA      NA      NA      NA
##                         MDS13   MDS14   MDS15   MDS16   MDS17   MDS18   MDS19
## Eigenvalue            0.03460 0.03337 0.03217 0.02941 0.02694 0.02620 0.02508
## Proportion Explained  0.02176 0.02099 0.02023 0.01850 0.01695 0.01648 0.01577
## Cumulative Proportion      NA      NA      NA      NA      NA      NA      NA
##                         MDS20   MDS21   MDS22   MDS23   MDS24   MDS25   MDS26
## Eigenvalue            0.02379 0.02246 0.02179 0.02106 0.02039 0.01971 0.01801
## Proportion Explained  0.01496 0.01413 0.01371 0.01325 0.01282 0.01240 0.01133
## Cumulative Proportion      NA      NA      NA      NA      NA      NA      NA
##                         MDS27    MDS28    MDS29    MDS30    MDS31    MDS32
## Eigenvalue            0.01620 0.013849 0.010275 0.009400 0.008019 0.007311
## Proportion Explained  0.01019 0.008711 0.006463 0.005913 0.005044 0.004599
## Cumulative Proportion      NA       NA       NA       NA       NA       NA
##                          MDS33     iMDS1
## Eigenvalue            0.006389 -0.004652
## Proportion Explained  0.004019  0.002926
## Cumulative Proportion       NA        NA
## 
## Accumulated constrained eigenvalues
## Importance of components:
##                        dbRDA1  dbRDA2
## Eigenvalue            0.05871 0.04871
## Proportion Explained  0.54653 0.45347
## Cumulative Proportion 0.54653 1.00000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  2.750505 
## 
## 
## Site scores (weighted sums of species scores)
## 
##           dbRDA1   dbRDA2     MDS1     MDS2       MDS3      MDS4
## SITE 1  -0.59427 -0.16451 -0.14771 -0.33607 -0.0859223 -0.212994
## SITE 11 -0.45337 -0.46011 -0.23978  0.15386  0.1261509  0.237370
## SITE 12 -0.54652 -0.75642 -0.44877  0.23840 -0.2252134 -0.115066
## SITE 13 -0.12157 -0.51390  0.05860 -0.62753 -0.0424431 -0.256505
## SITE 14 -0.31417 -0.28591 -0.04140 -0.38730  0.0396086 -0.203654
## SITE 15 -0.30246 -0.79945 -0.03439 -0.11317 -0.1309347 -0.211263
## SITE 16 -0.45766  0.40537 -0.02829 -0.07163  0.0308933  0.125735
## SITE 17 -0.50007 -0.43325 -0.06301  0.73800  0.4760396  0.361411
## SITE 18 -0.23591  1.09709  0.61113  0.80322 -0.2159329  0.452248
## SITE 19  0.55697  2.59188  1.51522  0.72623 -0.4623505  0.265834
## SITE 2   0.07451  0.03855 -0.38949  0.04156 -0.0701938  0.374170
## SITE 20  1.27192  0.05240  0.05244 -0.18703  0.8794177  0.804141
## SITE 21 -0.64222 -0.20393 -0.40958  0.09451 -0.0002606  0.117485
## SITE 22 -0.69477 -0.39862 -0.55866 -0.39024  0.5266838 -0.236789
## SITE 23 -0.23341 -0.07011  0.19956  0.31332 -0.3893083  0.086348
## SITE 24  0.41401 -1.20918 -0.35479 -0.23012  0.0081164  0.193247
## SITE 26 -0.53218  0.71759 -0.02716  0.45515 -0.0183075  0.359351
## SITE 28 -0.61852  1.00385 -0.02689  0.41929  0.0029756  0.430581
## SITE 29  0.47626  0.43365 -0.12521  0.12478 -0.5950115  0.342803
## SITE 3  -0.54249 -0.73219 -0.44976  0.05668 -0.1304354  0.046421
## SITE 30 -0.88529 -0.46309 -0.61315  0.33794 -0.2093706  0.150316
## SITE 31  0.15702 -0.98675  0.23967 -0.28872 -0.1116307 -0.524180
## SITE 32  0.06463  0.23715  0.51475 -0.29819  0.0294093 -0.097103
## SITE 33 -0.24919 -0.97106 -0.16670 -0.25156 -0.1300399 -0.272378
## SITE 34 -0.53685  0.34384 -0.58079  0.14271 -0.2847025  0.213550
## SITE 35  0.57845 -0.61546  0.17558  0.01290  0.5881978  0.538055
## SITE 36  1.13941 -1.04701 -0.16168 -0.55971 -0.5621101 -0.496346
## SITE 38  0.11771  0.53367  0.16337  0.22829 -0.3411253 -0.095943
## SITE 39 -0.71027 -0.83212 -0.87000 -0.10688 -0.2047559 -0.030265
## SITE 4   0.14585  0.39575 -0.22209 -0.06629 -0.2564930  0.011573
## SITE 40 -0.18800  0.14387  0.46789 -0.34817 -0.1578324 -0.142457
## SITE 41  0.07461  0.82427  1.00639 -0.38424  0.1709964 -0.674233
## SITE 42  2.33398 -0.53599  0.44914 -1.03948  1.4092795  0.827914
## SITE 5   0.27292  0.89196  0.20209  0.11649 -0.4067229  0.125238
## SITE 6   0.45087 -0.25704  0.32728 -0.67907 -0.2681291 -0.876440
## SITE 7   0.87161  1.21367 -0.16936  1.32868  1.3928171 -1.620007
## SITE 9   0.35845  0.81153  0.14558  0.03335 -0.3813597  0.001832
## 
## 
## Site constraints (linear combinations of constraining variables)
## 
##           dbRDA1    dbRDA2     MDS1     MDS2       MDS3      MDS4
## SITE 1  -0.53436  0.244983 -0.14771 -0.33607 -0.0859223 -0.212994
## SITE 11 -0.28142 -0.386202 -0.23978  0.15386  0.1261509  0.237370
## SITE 12  0.09204 -0.878470 -0.44877  0.23840 -0.2252134 -0.115066
## SITE 13 -0.48524  0.195319  0.05860 -0.62753 -0.0424431 -0.256505
## SITE 14 -0.54659  0.190041 -0.04140 -0.38730  0.0396086 -0.203654
## SITE 15 -0.16703 -0.628956 -0.03439 -0.11317 -0.1309347 -0.211263
## SITE 16 -0.68651  0.553341 -0.02829 -0.07163  0.0308933  0.125735
## SITE 17 -0.24577 -0.530833 -0.06301  0.73800  0.4760396  0.361411
## SITE 18 -0.28345 -0.436057  0.61113  0.80322 -0.2159329  0.452248
## SITE 19  0.39319  0.351596  1.51522  0.72623 -0.4623505  0.265834
## SITE 2   0.41166  0.294943 -0.38949  0.04156 -0.0701938  0.374170
## SITE 20  0.62134  0.131756  0.05244 -0.18703  0.8794177  0.804141
## SITE 21 -0.53883  0.168679 -0.40958  0.09451 -0.0002606  0.117485
## SITE 22 -0.62229  0.383342 -0.55866 -0.39024  0.5266838 -0.236789
## SITE 23 -0.01883 -0.585259  0.19956  0.31332 -0.3893083  0.086348
## SITE 24  0.51885 -0.483973 -0.35479 -0.23012  0.0081164  0.193247
## SITE 26 -0.53040  0.217589 -0.02716  0.45515 -0.0183075  0.359351
## SITE 28 -0.69213  0.575274 -0.02689  0.41929  0.0029756  0.430581
## SITE 29  0.48333  0.208739 -0.12521  0.12478 -0.5950115  0.342803
## SITE 3  -0.11176 -0.484151 -0.44976  0.05668 -0.1304354  0.046421
## SITE 30 -0.20601 -0.234679 -0.61315  0.33794 -0.2093706  0.150316
## SITE 31  0.02640 -0.958494  0.23967 -0.28872 -0.1116307 -0.524180
## SITE 32 -0.45895 -0.021464  0.51475 -0.29819  0.0294093 -0.097103
## SITE 33 -0.18938 -0.611963 -0.16670 -0.25156 -0.1300399 -0.272378
## SITE 34  0.18801  0.867602 -0.58079  0.14271 -0.2847025  0.213550
## SITE 35  0.20501 -0.526901  0.17558  0.01290  0.5881978  0.538055
## SITE 36  0.82790 -0.414978 -0.16168 -0.55971 -0.5621101 -0.496346
## SITE 38  0.53363  0.097754  0.16337  0.22829 -0.3411253 -0.095943
## SITE 39  0.40037  0.428319 -0.87000 -0.10688 -0.2047559 -0.030265
## SITE 4   0.32911  0.494131 -0.22209 -0.06629 -0.2564930  0.011573
## SITE 40 -0.44547  0.009707  0.46789 -0.34817 -0.1578324 -0.142457
## SITE 41 -0.51320  0.184444  1.00639 -0.38424  0.1709964 -0.674233
## SITE 42  0.57268 -0.023441  0.44914 -1.03948  1.4092795  0.827914
## SITE 5   0.47441  0.462943  0.20209  0.11649 -0.4067229  0.125238
## SITE 6   0.51953  0.314869  0.32728 -0.67907 -0.2681291 -0.876440
## SITE 7   0.51653  0.318764 -0.16936  1.32868  1.3928171 -1.620007
## SITE 9   0.44365  0.511687  0.14558  0.03335 -0.3813597  0.001832
## 
## 
## Biplot scores for constraining variables
## 
##         dbRDA1 dbRDA2 MDS1 MDS2 MDS3 MDS4
## netcen -0.2800 0.9600    0    0    0    0
## updist -0.9905 0.1372    0    0    0    0
# Generate a plot of the RDA
rownames(spe.rda$CCA$biplot) = c("network centrality","upstream distance")
ordiplot(spe.rda, choices=c(1,2), type="n", xlab="dbRDA 1 (54.65 %)", ylab="dbRDA 2 (45.35 %)")
points(spe.rda, "sites", pch=21, col="black", bg="white")
text(spe.rda, "cn", col="grey", cex=0.9, scaling=1, pos=4)
## Warning in arrows(0, 0, pts[, 1], pts[, 2], length = head.arrow, ...): "pos" is
## not a graphical parameter
## Warning in strwidth(labels, ...): "pos" is not a graphical parameter
## Warning in strheight(labels, ...): "pos" is not a graphical parameter

# Save the plot to an object
figure4b = recordPlot()

8.2.3 Variation partitioning

#Variation partitioning
spe.varpart1 <- varpart(meandist_bray, env_select[,1:8], env_select[,9:10])
plot(spe.varpart1,digits=2)

spe.varpart1
## 
## Partition of squared Unknown user-supplied distance in dbRDA 
## 
## Call: varpart(Y = meandist_bray, X = env_select[, 1:8], env_select[,
## 9:10])
## 
## Explanatory tables:
## X1:  env_select[, 1:8]
## X2:  env_select[, 9:10] 
## 
## No. of explanatory tables: 2 
## Total variation (SS): 1.5898 
## No. of observations: 37 
## 
## Partition table:
##                      Df R.squared Adj.R.squared Testable
## [a+c] = X1            8   0.28328       0.07850     TRUE
## [b+c] = X2            2   0.06757       0.01272     TRUE
## [a+b+c] = X1+X2      10   0.34269       0.08988     TRUE
## Individual fractions                                    
## [a] = X1|X2           8                 0.07716     TRUE
## [b] = X2|X1           2                 0.01138     TRUE
## [c]                   0                 0.00134    FALSE
## [d] = Residuals                         0.91012    FALSE
## ---
## Use function 'dbrda' to test significance of fractions of interest
anova.cca(dbrda(meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander + Condition(netcen + updist),
                data=env_select), step=1000)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander + Condition(netcen + updist), data = env_select)
##          Df SumOfSqs      F Pr(>F)  
## Model     8   0.4374 1.3603  0.013 *
## Residual 26   1.0450                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova.cca(dbrda(meandist_bray ~ netcen + updist+
                  Condition(T_av + O2_sat_av + Con_av + COD_av 
                     + NH4._av + Nt_av + pool_riffle + meander), data=env_select), step=1000)
## Permutation test for dbrda under reduced model
## Permutation: free
## Number of permutations: 999
## 
## Model: dbrda(formula = meandist_bray ~ netcen + updist + Condition(T_av + O2_sat_av + Con_av + COD_av + NH4._av + Nt_av + pool_riffle + meander), data = env_select)
##          Df SumOfSqs      F Pr(>F)
## Model     2  0.09446 1.1751  0.176
## Residual 26  1.04500

8.3 Figure 4

library(gridExtra)
#pdf("Figure4.pdf", width = 10, height = 10)
#grid.arrange(figure4a, figure4b, figure4c, figure4d)
#dev.off()

library(ggpubr)
library(gridGraphics)
## Loading required package: grid
pdf("Figure4.pdf", width = 16, height = 10)
ggarrange(plotlist = list(figure4a, figure4c, figure4b, figure4d), nrow = 2, ncol=2)
dev.off()
## png 
##   2

9. Structural equation model

9.1 PCA of environmental variables

# Conduct a PCA on the six physico-chemical variables from VMM
PCA = prcomp(environment2[, c(2:7)], center = TRUE, scale. = TRUE)
summary(PCA) # Summary statistics of the PCA
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6
## Standard deviation     1.5652 1.2707 0.9683 0.66323 0.62754 0.40523
## Proportion of Variance 0.4083 0.2691 0.1563 0.07331 0.06563 0.02737
## Cumulative Proportion  0.4083 0.6774 0.8337 0.90700 0.97263 1.00000
biplot(PCA) # Graphical representation of the PCA

var <- get_pca_var(PCA)
var$contrib # Show contribution of variables to PCs
##               Dim.1      Dim.2      Dim.3       Dim.4      Dim.5      Dim.6
## T_av       4.781229  3.2285483 86.3186845  1.24296167  3.7819969  0.6465796
## O2_sat_av 21.325906 12.6315198  4.5711032 37.97350069 10.9404721 12.5574980
## Con_av     1.259548 46.2188894  1.3713634 23.59374367 26.5347837  1.0216720
## COD_av    20.367872 15.5236356  1.6185615  9.07873660 46.8437587  6.5674354
## NH4._av   35.988850  0.3815385  0.8826774  0.05274591  0.3244479 62.3697404
## Nt_av     16.276595 22.0158684  5.2376100 28.05831146 11.5745407 16.8370745

9.2 Structural equation modelling

pc1 = PCA$x[,1]

data_p = data
levels(data_p$site) = as.factor(environment2$Site)
environment3 = cbind(environment2, pc1)
m = merge(x = environment3, y = data_p, by.x = "Site", by.y = "site", all.y = TRUE)
m2 = m[complete.cases(m[, c("SMI", "PI", "PI_ecto", "PI_endo")]), ] 
model_PI = psem(
  lme(PI ~ pc1 + SMI, random = ~1|Site, data=m2),
  lme(SMI ~ pc1, random = ~1|Site, data=m2)
  )
plot(model_PI)
summary(model_PI)
## 
## Structural Equation Model of model_PI 
## 
## Call:
##   PI ~ pc1 + SMI
##   SMI ~ pc1
## 
##     AIC      BIC
##  18.000   57.079
## 
## ---
## Tests of directed separation:
## 
##  No independence claims present. Tests of directed separation not possible.
## 
## Global goodness-of-fit:
## 
##   Fisher's C = 0 with P-value = 1 and on 0 degrees of freedom
## 
## ---
## Coefficients:
## 
##   Response Predictor Estimate Std.Error  DF Crit.Value P.Value Std.Estimate  
##         PI       pc1   0.2698    0.1154  35     2.3376  0.0253       0.1687 *
##         PI       SMI   0.1889    0.3794 530     0.4980  0.6187       0.0199  
##        SMI       pc1   0.0040    0.0088  35     0.4590  0.6491       0.0239  
## 
##   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05
## 
## ---
## Individual R-squared:
## 
##   Response method Marginal Conditional
##         PI   none     0.03        0.17
##        SMI   none     0.00        0.04
model_PI_ecto = psem(
  lme(PI_ecto ~ pc1 + SMI, random = ~1|Site, data=m2),
  lme(SMI ~ pc1, random = ~1|Site, data=m2)
  )
plot(model_PI_ecto)
summary(model_PI_ecto)
## 
## Structural Equation Model of model_PI_ecto 
## 
## Call:
##   PI_ecto ~ pc1 + SMI
##   SMI ~ pc1
## 
##     AIC      BIC
##  18.000   57.079
## 
## ---
## Tests of directed separation:
## 
##  No independence claims present. Tests of directed separation not possible.
## 
## Global goodness-of-fit:
## 
##   Fisher's C = 0 with P-value = 1 and on 0 degrees of freedom
## 
## ---
## Coefficients:
## 
##   Response Predictor Estimate Std.Error  DF Crit.Value P.Value Std.Estimate    
##    PI_ecto       pc1   0.3092    0.0851  35     3.6342  0.0009       0.2812 ***
##    PI_ecto       SMI   0.1335    0.2525 530     0.5286  0.5973       0.0205    
##        SMI       pc1   0.0040    0.0088  35     0.4590  0.6491       0.0239    
## 
##   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05
## 
## ---
## Individual R-squared:
## 
##   Response method Marginal Conditional
##    PI_ecto   none     0.08        0.25
##        SMI   none     0.00        0.04
model_PI_endo = psem(
  lme(PI_endo ~ pc1 + SMI, random = ~1|Site, data=m2),
  lme(SMI ~ pc1, random = ~1|Site, data=m2)
  )
plot(model_PI_endo)
summary(model_PI_endo)
## 
## Structural Equation Model of model_PI_endo 
## 
## Call:
##   PI_endo ~ pc1 + SMI
##   SMI ~ pc1
## 
##     AIC      BIC
##  18.000   57.079
## 
## ---
## Tests of directed separation:
## 
##  No independence claims present. Tests of directed separation not possible.
## 
## Global goodness-of-fit:
## 
##   Fisher's C = 0 with P-value = 1 and on 0 degrees of freedom
## 
## ---
## Coefficients:
## 
##   Response Predictor Estimate Std.Error  DF Crit.Value P.Value Std.Estimate 
##    PI_endo       pc1  -0.0538    0.0792  35    -0.6792  0.5015      -0.0454 
##    PI_endo       SMI   0.1217    0.2894 530     0.4205  0.6743       0.0173 
##        SMI       pc1   0.0040    0.0088  35     0.4590  0.6491       0.0239 
## 
##   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05
## 
## ---
## Individual R-squared:
## 
##   Response method Marginal Conditional
##    PI_endo   none        0        0.11
##        SMI   none        0        0.04